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A COMPARISON OF AN ANALYTICAL AND NUMERICAL 
SOLUTION OF THE RELATIVE MOTION 
OF TWO CLOSE SATE LUTES OF 
AN OBLATE PLANET 

Richard W* Barbieri 
Joel J, Mashbaum 


ABSTRACT 

This report presents the results of a study undertaken to compare an ana- 
lytical theory of relative motioUj developed by R. Barbieri (Ref. 1) , with a numer- 
ical integration of the equations of planar relative motion. 

The equations of relative motion were derived from a Lagrangian formula- 
tion and is presented in detail in [1]. The numerical integration has been carried 
out with a program that integrates a system of ordinary differential equations 
which may or may not be coupled. On the other hand^ the analytical solution has 
been constructed successively in terms of powers of the eccentricity of the ref- 
erence satellite and powers of a small parameter. Since for the application to 
problems of orbiting long baseline interferometry the eccentricity is small 
(< .01), the expansion of the solution in powers of eccentricity is carried to the 
first power only. 

The results indicate, as ejq>ected, good agreement, using small initial ve- 
locities, over several days for a range of eccentricities. For example, at a 
semi -major axis of 17,000 km, and eccentricity of .01 and initial velocity com- 
ponents of .Olkm/secthe analytical development yields a baseline distance which 
differs from the distance as calculated numerically by no more than 5% over the 
first 3 1/2 days. On the other hand when the semi-major is increased to 20,000 
km, the eccentricity decreased to .0001 the baseline distance as calculated by 
both programs differ by no more than 5% over tiie first 6 days. Comparisons 
are presented to show the effect of including the oblateness of the central body. 
Results are documented to indicate the dependence of the motion of one satellite 
about the other upon a variety of initial conditions. 
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GLOSSARY OF SYMBOLS 

radius vector of the M satellite with respect to the center of 
mass of the central body 

radius vector of the O* satellite with respect to the center of 
mass of the central body; 

position vector of M with respect to O®; has components 

{x» y, 2} 

the semi-major axis of the orbit of O^ 
the eccentricity of the orbit of 
the inclination of the orbit of 

the longitude of the ascending node of the orbit of O' 
the argument of perigee of O' 

the secular rate of change of the argument of perigee 
the true anomaly of O’ 
the sum (w + v) 

the angular velocity of O’; has components {c<j » oj #0; } 

K y z 

the velocity components of O’ in the rotating coordinate system 

the components of force, in the rotating coordinate system, acting 
on the M satellite. 

the gravitational constant of the central body 

the equatorial radius of the central body 

the first non zero coefficient in the spherical harmonic expan- 
sion of the gravity field of the central body 

the potential energy of the M satellite 








A COMPARISON OF AN ANALYTICAL AND NUMERICAL 
SOLUTION OF THE RELATIVE MOTION 
OF TWO CLOSE SATELLITES OF 
AN OBLATE PLANET 


INTRODUCTION 

In Reference [1] a solution of the equations of planar relative motion was 
constructed taking into account the oblateness of the central body. Since the 
publication of those results an effort has been made to determine how the ana- 
lytical solution compares to the numerical integration of the equations of planar 
relative motion. This report represents a culmination of that effort. 

A study of the relative motion of two close spacecraft orbiting an oblate 
planet has some rather interesting applications. For example the concept of an 
orbiting long baseline interferometer has some intriguing implications. Not 
only is it possible to get data unperturbed by the troposphere or ionosphere 
from various radio sources but it is conceivable that such a system, yielding 
extremely accurate positions of radio sources, can be used to provide an inde- 
pendent estimation of the structure of the troposphere and ionosphere and con- 
sequently an independent verification of refraction effects. 

The analytical representation was constructed by expanding the equations 
of motion into powers of eccentricity and a small parameter, e , which is related 
to the second zonal harmonic coefficient . The equations of motion resulting 
from this operation are then solved by successive approximations. The first 
step in this successive approximation technique yields an almost periodic solu- 
tion, independent of e and eccentricity, whose frequency depends upon the modi- 
fied mean motion. The second step introduces the dependence of the motion 
upon the eccentricity of the reference satellite and also yields an almost peri- 
odic solution with three frequencies appearing. The last step of the procedure 
accounts for the obiateness of the central body and yields a solution which de- 
pends upon four frequencies, Uiree of which are dependent upon the secular rate 
of increase of the argument of perigee. Furthermore the last step of the pro- 
cedure also Introduces a purely secular term; mixed secular terms have been 
eliminated by making the angular velocity of the reference! satellite dependent 
upon the small parameter. 

The routine FNOL2 uses the fourth order Runge-Kutta and/or Adams- 
Moulton methods of solution of ordinary differential equations which may or 
may not be coupled.' As currently modified, the program uses a double precision 
arithmetic throughout. Furthermore the truncation error can be held within 
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Input bounds by giving the user an option of automatically varying the step size. 
More detailed information is available in [2] , 

A description of the subroutines used in the numerical and analytical pro- 
grams is presented in Appendix n and Appendix HI respectively, A listing of 
the subroutines is given in Appendix IV and Appendix V. 


THE EQUATIONS OF RELATIVE MOTION 


The equations of relative motion which are of main concern in this report 
are derivable from the following system of coupled differential equations! 


X + - 2 y Ctjj, - y ” ^o' (z - x 

m 3x - 

y +yg' + 2 (xoj^-z + X - z + Xq/ 


■ 2 2 f ® ^ c 


z + Zq. + 2 y + y + y^* (z 


X ^ 


.i |“ = F, 

m 3 z 

See Figure 1 for the geometry. 
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A solution to this system of equations presents a formidable task since they 
are coupled. In Reference [l] the angular velocity component a ^is neglected 
thereby uncoupling the motion; this imposes a constraint on the rotation of the 
orbital plane of O* about , In particular the motion governed by the differential 
equation I cos u + ^ siii u sin l is being neglected. Making the substitutions 


U = - 




( 4 ) 


■3 X -u Xnf 

- 5 “ h a 

x5, 


- 3 s in^ L s in* u) 


( 5 ) 


^ 

'‘o' 


( 6 ) 


Xo' -«o' 


xh 


JL i (X*« &: ) = (F , 

Xjj. d t 0 r' ^ y'o 


( 8 ) 


reduces the equations of planar relative motion to the form 


X - 2yct;^-ya)^-xa)* + Jt. * 
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These represent the equations which have been progframmed for numerical so- 
lution by FNOL2, 

In order to present an analytical solution of these equations certain trans- 
formations and representations must be used so as to have the time appearing 
explicitly. To begin, Is transformed to 


5/ 


6 ' 


5 


~ 2 — — —I e s in V 

a (1 - p2) 


and then the true anomaly of the reference satellite Is introduced via the trans- 
formation 


a (1 - e^> 

1 e cos V 


( 11 ) 


Finally the true anomaly is expressed as a function of the mean anomaly through 
the representation 


) exp (i > v) 


exp i M (a - 2 >) exp i (y + 1) 


( 12 ) 


w. 


(a + 2 y) exp i ( v - 1) M 


where 


tt, y = 0, + 1, ± 2, ... 

Making the necessary substitutions leads to the following system of equations 
which are more amenable to analytical solution 


u 'll n, . 2 e sin V 9 

x-2yaj_+2yct); x 


e cos V 


1 + e cos V “ 1 + e cos v 


~ ^ X (1 - (1 + e cos v) (1 « 3 sin^ t sin^ u) 


( 13 ) 
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+ 2 X - 2 X Ci-2 


e sin V 
1 -r cos V 


3 e cos V 

y uji 

* 1 t e cos V 



(14) 
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These equatiom;, as has been mentioned, are solved using the method of 
successive approximations starting with a representation of the solution in the 
form 


X (t) ^ (t) + e Xqj (t) + e XjQ (t) 

y (t) = (t) + e (t) f € (t) 

Substituting these two expressions into (13) and (14), using (12) to get the mean 
anomaly and expressing in the form 

cli =n+Ocos 1 , + 2necosv + £ii (17) 

* 2 

yields the system of equations which govern the first order behavior of the 
satellite I 


(15) 

(16) 


\o - 2 " Voo = 0 

*01 ■ 2 ff y„, - 4 n cos n t + 2 n* y„„ sin n t 
- XjjQ cost! t = 0 



- 2 n 
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sin^ t cos (2 + 2 n) t 
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^00 + 2 n ^ 0 


Vqj + 2 n Xqj + 4 n cos n t - 2 n ^ s in ti t 


- y cos n t - 0 


/R \2_ 

^10 - 2 n Xj^ T 2 x^o - I— ) »2 sin2 , cos (2 w + 2 n) t. 

\ 3 / 


where the expressions M = nt and w = t have been used. The solution of the 
system (18) - (19) is presented in [1] ; the results are that 


Xqq ( t) Uj sin 2 n t + a cos 2 n t 


(tl = cos 2 n t - a- sin 2 n t j- 


sin n t + cos n t + /iJJ sin (2 n + n) t 
+ /3JJ sin (2 n - n) t + cos (2 n + n) t 
+ cos (2 n - n) t 


y®! (t) = sin n t + ^os n t + sin (2 n + n) t 
+ sin (2 n - n) t + cos (2 n + n) t 


+ cos (2 n - n) t 


’^lo 2 (Wq -u n) t + /3JJ cos 2 (w^, + n) t 

+ sin 2 (Wq + n + n) t + /5JJ sin 2 (w^ + n - n) t 
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+ 311 cos 2 (Wq + n + n) t + /3jo cos 2 (w^ + n - n) t 


+ /3^o. 

^ ■ 7x* 


yio = >^1; sin 2 n t + ^10 cos 2 n t + sin 2 (w„ + n) t 
+ cos 2 (w„ + n) t + sin 2 (w„ + n + n) t 

+ sin 2 (w, + n - n) t + ,3“ cos 2 (w„ + n T n) t 

+ cos 2 (w„ + n - n) t + ,S>; t. 


(25) 


■“ “i. 3 = »• 1. 2. 3, {fi%. ^»J> With j = X. 2 6 and 

^Pjx» /3jy>with ] = 1, 2 7 are given explicitly in the Appendix I. 

Turning attention to the behavior of the motion of the reference satellite, 
it will be noticed that in the numerical integration program this behavior is 
specified by a modification of the planetary equation of Lagrange to allow for 

veiy sm^ eccentricities. To be specific, the Lagrange equations, before any 
modifications, take the form ^ 


. — {S e sin v + T (1 + e cos v)} 

n N 1 ^ 
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'll - 
fTa 


S 


cos V 



where S is the force due to the first zonal harmonic directed along the radius 
vector to O*, T is the force due to the first zonal harmonic 90® ahead of S and 
located in the orbit plane of 0», W completes the right handed coordinate system. 

Because of the presence of the eccentricity in the denominator of the eo and 

^ motion for very smaU eccentricities becomes quite 

difftcult. To overcome this situtation a transformation of variables is made; in 
particular let, 

e cos w 
n» = e sin w 

and replace the equation for v by the equation for u where u = v + w. These 

tr^formations lead to the following system of equations describing the motion 
of the reference satellite: "iwwon 




~ m cos u) 


1 - 3 sin® i sin® u 
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> V 


in terms of the mean anomaly, M, is presented. The Fourier development of 
<r/ a) ^ exp i r v is given as 



exp i y V ~ 



exp i k M 


(28) 


where the development shows that the coefficient independent of M (and hence 
independent of time explicitly) is given by 


( „ 1)7 r 2) (u 3). ♦ y -j- 1) 


( 


y - a - I 7-a 


1 - y, e 




(29) 


that is, is the constant coefficient in the Fourier development; the depend- 
ence on the eccentricity is explicitly shown in the expansion - 

Carrying out these operations in the d equation of (26) leads to an equation 
tor the secular variation of w, namely w^, given by 


'*'0 =1*' Jo — C5 cos2 t -1), 


R? 


(30) 


The remaining set of equations governing the motion of the reference 
are given by 


satellite 


a--3Mj. 


R? 


*0" 


1 

n >fl ^ e* 


*^(1 - 3 sin* 


I sin* u) e sin (u - 





sill? V sin s 2 


“} 



• 3 _ nI 1 ^ 02 

— — (l-3sin2 L sin2 u) sin (u - w^) 




“V (" - %> + 2 cos (u - «- ) + e] sin* i sin 2 u 


I = 


3 /^J 2 ^e si n 2 l sin 2 u 

^ n a2 41 _ q2 


i 

*o' 


(31) 



• tan u 
sin t 


u = 


n Br 'll — e® 

X?f 


- n cos 


where 



a (1 - e^) 

+ e cos (u - Wq) 


THE RESULTS 


nti. ti. ** ** compared are the reaidta generated by (20) thru (25). 

with the results of the numerical integration of (^ and (10). The purpora of 

s comp on is to determine the Interval of time over which the anatotioal 
solution is valW and to illustrate the dependence of the analytical develf^uneto 
•m smaU eccentricities. Furthermore it wlU be shown just how critical the 
initial relative velocity is to the motione 

Compart^ h^ been obtained for about 25 cases where the semi-major 
^ ^ 12.000 tan and 20,000 tan, toe eccmitriclly varies be- 

tween .0001 and.Oland toe initial relative velocity varies between 1 m/sec to 


10 m/sece All cases have the inciination of the reference satellite initially at 
10°, the longitude of the asceviding riode at 20°, the argument of perigee at 0?0, 
and the true anomaly initially at 45®o 

Tables I and II were generated for a semi -major axis of 17,000 km and 
eccentricity of »01; the only parameter which was varied was the initial velocity. 

It is noticed that an increase of one order of magnitude in the velocity increases 
the baseline magnitude by a factor of roughly 10; coupled with this increase is a 
significant difference in the orientation of the baseline. 

Table I'l differs from Table I in that only the eccentricity is different. As 
expected, the difference in the magnitude and orientation of the baseline is 
negligible. 

Tables IV and V were generated with a semi-major axis of 20,000 km, and 
a nearly circular orbit; once again only the initial velocity was changed. It is 
noticed immediately that the baseline undergoes slightly larger excursions when 
compared to the cases where a = 17,000 km; furthermore the orientation of the 
baseline is much different due to the presence of n andn* in the equations (20) 
thru (25), The agreement between the analytical and numerical is much better 
than when the semi- major aads is 17,000 km; in fact not only is the agreement 
better but it lasts longer. For example in Table I it is noticed that after about 
5 days the agreement between the x components begins to break down. 

Table VI shows what can hapi>en with a slight variation in the initial velocity. 

Finally the graphs I through IV are Included to show the accumulative 
effect of the oblateness of the central body, in these cases the earth. What is 
shown is the variation of the difference in baseline length versus time, and 
are the baseline lengths as computed in the analytical program and numerical 
program respectively. The insensitivity of the motion as a function of eccentricity 
is clearly shown as is the strong dependence of the motion on the semi -major 
axis of the reference satellite. 


COMMENTS, CONCLUSIONS AND SUMMARY 

What has been presented are some results generated by comparing an 
aiialytical development of the equations of relative motion to a numerical solu- 
tion of these equations. The results clearly indicate the motion is strongly 
dependent upon the semi-major axis of the reference satellite and on the Initial 
relative velocity. For example the following summary shows the dependence 
on the aforementioned parameters and indicates the percent difference between 
the analytical value of the baseline distance and the numerical value; also shown 
is the interval of time the two values agree to within 5 percent. 
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a (km) G 


km/s km/s 


Time 

Interval 


% Difference of 
Baseline Length 


17,000 

.01 

.01 

.01 

3 1/2 Days 

< 5 

17,000 

.01 

.001 

.001 

3 1/2 Days 

< 5 

17,000 

.001 

.01 

,01 

3 l/2 Days 

<5 

20,000 

,0001 

.01 

.01 

6 Days 

<5 

20,000 

.0001 

.001 

.001 

6 Days 

<5 


Qualitative and quantitative Information has been presented to indicate the 
accumulative effects of oblateness upon the relative motion; it is clear from the 
cases shown that a variation of the eccentricity from 0,0001 to 0.01 has little 
effect on the relative motion over a period of 7 days. For example the following 
table presents the difference between the numerical and analytical values for the 
baseline length at the 7th day; the initial components of velocity are ,001 km/s. 


a (km) 

e 

(Ra - R„) km 


17,000 

.0001 

- .1538 

off 

17,000 

.0001 

+ .167 

on 

17,000 

.001 

- .1680 

off 

17,000 

,001 

+ .151 

on 

20,000 

.0001 

- .0114 

off 

20,000 

.0001 

- .312 

on 

20,000 

.001 

+ .0063 

off 

20,000 

.001 

- .2932 

on 


A time history over 7 days of the difference - Rj^ is shown in the griq)hs. 

Agreement of the analytical solution with the numerical is much better at 
larger semi-major axes since the periodic behavior of Q and w is less pronounced 
at the semi-major axes of 20,000 km or more. 

It must be mentioned that the analytical development was carried out under 
the assumption that n and h are constant. Being rigorous, this assumption is 
not quite valid; the following two tables show how n, n and the coefficients of 
integration vary over an interval of seven days. 
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a = 20,000 km, e = ,001 Jj 0 
,2232143 X 10"^ /sec < n ^ .2232177 x lo”^ /sec 
,2231463 X 10"^ /sec < n < ,2232145 x 10"^ /sec 
,143879 X 10"^ km/sec ^ .144600 x 10“^ km/sec 

2.233363 km ^ <2,238736 km 

- 2.237675 km £ ^ - 2,238756 km 

- 1.236651 km ^ ^ - 1,237939 km 


a = 17,000 e = .001 / 0 

.2848343 x 10“^ /sec < n ^ .2848404 x 10"^ /sec 
.2847141 X lO’^/sec i n ^ ,2848344 x lO'^/sec 
.1569019 X 10*2 km/sec ^ ^ .1569368 x lO"^ km/sec 

1,753537 km ttj ^ 1.754718 km 
“ 1,753565 km ^ ^ “• 1,754745 km 

A 

- 0,7525190 km^ <^3 ^ ' 0.7539899 km. 

Comparisons of. the use of w in the numerical program versus the use of 
Wo in the analytical program indicate that for a semi-major axis of 17,000 km 
the analytical value of u = v + Wo differs from the numerical value of u « v + w 
by about one to three degrees. The parameter was used in the analytical 
development to avoid consideration of differential equations with almost periodic 
coei^icients which arise when the oblateness Is included. The agreement shown 
here between the analytical and numerical programs indicates that the former 
would serve as an excellent tool for very early orbit determination of the rela- 
tive motion of two close satellites. For the studies alluded to in the intro- 
duction, namely, the study of spectra from radio sources and, the study of 
continental drift, a very precise orbit determination program is required. This 
means includii^ the perturbation effects chie to higher order harmonics and the 
presence of the sun and moon. Naturally then, an analytic solution of the form 
presented in [1] would be quite cumbersome and time consumii^ to carry out. 
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TABLE I 


a = 17,000 Xq = 1.0 km x^, = ,01 km/s Jj ^ ® 
e = ,01 ^0 ~ km/s 


'' ~ ' " ' ~"‘‘ ' n 

Time (hrs.) 

- -J 

X (km) 

! y (km) 

R (km) 


Numerical 

Analytical 

Numerical 

Analsrtical 

Numerical 

12 

-5.3791 

-5.6772 

-10.222 

-9.9053 

11.551 

11.415 

24 

-5.1736 

-6.0541 

-23.046 

-22,872 

23.619 

23,660 

36 

1.4709 

.1904 

-34.015 

-34.464 

34.046 

34.465 

48 

12.812 

11.552 

-40.184 

-41.484 

42.178 

43,062 

60 

25.854 

25.078 

-39.748 

-41.897 

47.417 

48.829 

72 

36.934 

37.148 

-32.622 

-35.429 

49,278 

51.335 

84 

42.745 

44.465 

-20.745 

-23.675 

47.513 

50.375 

96 

41.555 

44.957 

-7.6492 

-9,6819 

42.253 

45,988 

108 

33.913 

38.346 

2.8116 

2.8599 

34.030 

38,452 

120 

22.307 

26.221 

7.8796 

10.594 

23.658 

28.281 

132 

10.026 

11.633 

6.6392 

11.377 

12.025 

16.272 

144 

.04627 

-1.6942 

-.15822 

4.8340 

.1648 

5.1223 

156 

-5.5229 

-10.306 

-10.720 

-7,5457 

12.059 

12.773 

168 

-5.5038 

-11.886 

-22.719 

-22,760 

23.376 

25,677 




■■■ 










a (hrs.) 


Numerical 

12 

.36240 

.3333 

24 

.38386 

.2960 

36 

1.0497 

.9201 

48 

2.1858 

2.0557 

60 

3.4922 


72 

4.6025 

4.6153 

84 

5.1859 


96 

5,0690 

5.4012 

108 

4.3066 I 

4,7449 

120 

3.1473 i 

3.5384 

132 

1.9200 

2.0857 

144 

.9222 1 

,7584 

156 1 

,3649 j 

-.0987 


168 .3660 -.2545 












a = 17,000 

B = .001 


TABLE m 


Xq = 1.0 km Xq 
y„ = 1.0 km 


.01 km/s 
,01 km/s 


Time (hrs.) | 

X (km) 

Analytical 

Numerical 

12 

-5.4330 

-5.7265 

24 

-5,3963 

-6.26/6 

36 

1.0243 

-.2816 

48 

12.203 

10.842 

60 

25.264 

24.313 

72 

36.598 

36.636 

84 

42.819 

44.512 

96 

41.988 

45.725 

108 

34.453 

39.748 

120 

22.672 

27.941 

132 

10.097 

13.236 

144 

-.0809 

- .5861 

156 

-5.6315 

-9.8678 

168 

-5.3791 

-12.038 1 


-10.196 

-23.018 

-34.106 

-40.543 

-40.459 

-33.621 

-21.800 

-8.4619 

+2.4078 

7.8078 

6.6514 

-.3015 

-11.119 

-23.316 


Numerical 


-9.8802 

-22.842 

-34.553 

-41.893 

-42.799 

-36.851 

-25.427 

-11.394 

1.6073 

10.105 

11.713 

5.7899 

-6.3433 

-21.695 


B (km) 


Analytical 

Numerical 

11.553 

11.419 

23.642 

23.686 

34.121 

34.555 

42.340 

43.273 

47.699 

49.223 

49.697 

51.963 

48.049 

51.263 

42.832 

47.123 

34.537 

39.781 

23,978 

29,712 

12,091 

17.674 

.3122 

5,8194 

12.464 

11.731 

23.929 

24.811 
















TABLE IV 


a - 20,000 km = 1.0 km 

e = .001 ~ 1.0 km 


Time (hrs.) 

! 

1 X (km) 

y 

Analytical 

Numerical 

Analytical 

12 

12.439 

12.660 

8.3138 

24 

25.827 

26.500 

10.169 

36 

38.758 

39.816 

6.2874 

48 

48.921 

50.074 

-2.6480 

60 

54.447 

55.232 

-15.092 

72 

54.228 

54.290 

-28.825 

84 

48.133 

47.331 

-41.266 

96 

37.261 

35.719 

-49.970 

108 

23.674 

21.564 

-53.150 

120 

10.110 

7.6823 

-50.249 

132 i 

-.8132 

-3.3527 

-41.913 

144 1 

-7.1295 

-9.3127 

-29.965 

156 

-7.9050 

-9.0567 

-16.673 

168 

-3,2543 

-2.5861 

-4.4052 


U 

yc 

^ = .01 km/s 

J ^0 
2 



[ R (km) 


Numerical 

Analytical 

Numerical 


8.4921 

14.962 

15,244 


10,286 

27.757 

28.427 


5.9627 

39.265 

40.260 


-3.6902 

48,992 

50.210 


-16.867 

56,500 

57.750 


-31.037 

61.413 

62.536 


-43.513 

63.399 

64.293 


-51.846 

62,333 

62.959 


-54,426 

58.184 

58.543 


-50.722 

51.257 

51,301 


-41,372 

41.921 

41.507 


-28.210 

30.801 

29,707 


-13.669 

18,452 

16.397 


-.6384 

5,4769 

2,6637 




















a = 20,000 km x^j 

6 = .0001 yjj 


X (km) 


Time (hrs.). 

Analytical 

Numerical 

! 

12 

2.1438 

2.1665 

24 , 

3.4827 

3.5514 

36 j 

4.7758 

4.8835 

48 

5.7921 

5.9102 

60 

6.3448 

6.4264 

72 

6.3228 

6.3330 

84 

5.7135 

5.6376 

96 

4.6260 

4.4772 

108 

3.2675 

3.0623 

120 

1.9108 

1.6748 

132 

.8187 

.5720 

144 

.1867 

t -.0231 

156 

,1095 

.0033 

168 

.5743 

.6510 


— '4 




<l',U 




TABLE V 


= loO km Xq = .001 km/s 

= 1.0 km Jq = .001 km/s *^2 ^ ® 


y (km) 

R (km) 

Analytical 

Numerical 

Analytical 

Numerical 

1.7308 

1.7477 

2.7553 

2.7835 

1.9163 

1.9260 

3.9751 

4.0401 

1.5276 

1.4920 

5.0142 

5.1064 

.6338 

.5255 

5.8266 

5,9335 

-.6109 

-.7937 

6.3742 

6.4753 

-1.9846 

-2.2119 

6.6269 

6.7082 

-3.2285 

-3.4610 

6.5626 

6.6152 

-4.0999 

-4.2955 

6.1814 

6.2046 

-4.4180 

-4.5549 

5,4950 

5.4887 

-4.1286 

-4.1858 

4.5494 

4.5085 

-3.2951 

-3.2521 

3.3952 

3,3020 

-2.1009 

-1.9373 

2.1091 

1.9374 

-.7718 

-.4845 

.7795 

.4845 

.4544 

.8169 

.7323 

1.0446 




' !• 











TABLE VI 


a = 20,000 km = 1.0 = .001 

e = .0001 Yq = 1.0 Yq = «01 


mm 

lll^llllllllllll^lllll^l^^^^^Q 

y (km) 

R (km) 

Analytical 

Numerical 

Analytical 

Numerical 

Anal 3 rtical 

Numerical 

12 

3.9966 

4.0336 

10.166 

10.359 

10,924 

11.116 

24 

10.520 

10.830 

17.208 

17.565 

20.169 

20.636 

36 

19.376 

20.092 

20,887 

21.170 

28.490 

29.187 

48 

28.982 

30,070 

20.542 

20.469 

35.524 

36,375 

60 

37.617 

38.828 

16.181 

15.534 

40.950 

41.820 

72 

43.672 

44.678 

8.5252 

7.3078 

44.497 

45.272 

84 

45.913 

46.456 

-1.0618 

-2.6941 

45.925 

46,534 

96 

43.848 

43.812 

-10.748 

-12.511 

45.146 

45.563 

108 

37.802 

37.180 

-18.614 

-20.308 

42.136 

42.365 

120 

29.032 

27.871 

-23.127 

-24.527 

37.118 

37.126 

132 

19.269 

17.586 

-23.462 

-24.359 

30.361 

30.044 

144 

10.414 

8.3780 

-19.737 

-19.811 

22.316 

21.510 

156 

4.0032 

1.9678 

-12.779 

-11.708 

13.392 

11.872 

168 

1.0361 

-.3454 

-3.9434 

-1.6383 

4.0772 

1.6743 




i 
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APPENDIX I 


The functions appearing (20) - (25) are 



X 

(0) 


y (0) 



D, . 

-Bi 

C 2 

^2 

y 

(0) 

-Di 

X (0) 



Cl 

*■ ^2 


A. 

* 

y 

(0) 

-C2 

, X (0) 


A. 

I>2 

-B. 

C 2 

5* 

■ 

X 

(0) 

- A. 

;y (0) 



Cl 

-A, 



where 


Ai 


1 n 4 n ~ n 

- e 


2 n 2 n 4 ^^2 _ ,,2 


B. 




3 
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n sin^ t 

. .i (5 


-f ^ 

16 \a 
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2 16 \a^ 

,n 4 
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16 (w^ 
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3 (Wq * n) 


+ 5 n ~] 
n + w^j + n ) 


PP^ING PAOB BUNK NOT PTIAfFn 


s i 


29 


n 8 n ^ - n n - 16 - < * 

2n(2n+n)(2n-n) \rf 



1 

( + n r FT) ^ 


('^0 



£ 

16 



2 



Sin 


2 I 


+ n 


3(wQ+n)n + 5n^ 

(Wq + n + FD^ (2 n + Wq + n) 



3(wQ+n^n-5n^ 

+ n - FF)^ [2 Fi - (Wq + n)] 


4 n - n ^ sin^ t ri^ - 2 (w^ + n)^ 
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A 
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= 2 n 


+ e 


2 n 2 
n 
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+ 
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n 
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n 

3 n - 8 n ^ 

2^ 
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n 
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n 

3 n - 8 n 
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The coefficients of Equation (23) are 
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The coefficients of Equation (25) are 
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NUMERICAL PROGRAM 


General Overview 


The numerical program Is comprised of a MAIN routine and six (6) sub- 
routines. The MAIN routine reads in the data (via 2 namelists) and writes out 
the initial conditions. Execution starts by integrating the equations of motion 
for the desired time period. Values for the orbital elements are written for 
each time point. Upon completion of the integration the values for position and 
velocity are plotted vs, time on the SD 4060, 
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MAIN Routine Description 


Purpose; To store initial values read in of the variables to be integrated. The 

variables not read in are calculated from input variables. The following 
variables are integrated; aje ,fi,x,y,z,x ,y, z. 

■i and m are obtained from e and w by; 

't- e cos w 
m = e sin w 

u is defined as u = v + w. 

The integrator (FNOL2) is then called. All values are calculated for each 
time point until the specified final time is reached. 



Subroutine DERIV Description 

Purpose; To compute values of the differentials to be used by the Integrator 
(FNOL2) 

Method; The systems of differential equations {a, igfl , I ^ m# u } and ( x, y, z» 
X, y, z> are calculated from Equations (27), (9) and (10). The inter- 
mediate variables used in the Integration of these equations are ob- 
tained from expressions for 

P " 

« 

OJ 

X 


a; 

X 

These expressions are to be found in the listing of the numerical program in 
Appendix 4. 





I 

I 


Subroutine Description 


h MAIN 

Purposes Tq read in control parameters and initial conditions. 

To write out control parameters and initial conditions. 

C alias DATK 
FNOL2 
PLOT 

Methods All input through use of namelists 

COMMON blocks useds INCOND 

PLOTS 

CONST 

FLAG 


Variables not in COMMONS 


FORTRAN Name 

Format 

Description 

A 

R*8 

semi -major axis 

E 

R*8 

eccentricity 

I 

R*8 

inclination 

J 

1*4 

no. of time points 

R (500) 

R*8 

position 

T (500) 

R*3 

time 

V 

R*8 

true anomaly 

w 

R*8 

argument of perigee 

X 

R*8 

X-coord, 

Y 

R*8 

Y- coord. 

Z 

R*8 

Z -coord. 

DR (500) 

R*8 

velocity 

NO 

1*4 

*NO» 

PL 

1*4 

control parameter to PLOT or 
not to 

PR 

P4 

control parameter to print 
Intermediate computed values 

TF 

R*8 

final time 

TI 

R*8 

initial time 

XO 

R*8 

distance between origins of 2 
coordinate systems (XO on fig.) 
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FORTRAN name 


Format 


Description 


FMU 

R*8 

- 3,986032 D5 

J20 

R*8 

Jjo = 1.0823D-3 

RSQ 

R*8 

r^ = 4.068098 D7 

YES 

1*4 

'YES’ 

DELM <20) 

R*8 

derivative array 

ELEM (20) 

R*8 

Integrated array 

STEP 

R*8 

integration interval 

XDOT 

R*8 

X-component of velocity 

YDOT 

R*8 

Y-comronent of velocity 

ZDOT 

R*8 

Z -component of velocity 

EPCXJH 

R*8 

epoch date in form: YYMMDD.D 

KPLOT 

L*4 

plot or not 

OMEGA 

R*8 

longitude of ascending node 

KOMEGA 

1*4 

parameter to set = 0, 

KPRINT 

1*4 

parameter to print intermediate 
values 

KUTMOL 

1*4 

type of Integration scheme to use 


41 


II. 


DERIV 


Purpose? To define variables to be integrated 
Called by; FNOL2 

Method? Variables to be integrated in equations in terms of predefined 
variables. 


COMMON blocks used? CONST 

FLAG 
ALFA 
A NS 


Variables not in COMMON? 


FORTRAN name 

Format 

A 

R*8 

C 

R*8 

E 

R*8 

I 

R*8 

N 

R*8 

P 

R*8 

T 

R*8 

U 

R*8 

V 

R*8 

X 

R*8 

Y 

R*8 

A1 

R*8 

A2 

R*8 

B1 

R*8 

B2 

R-8 

Cl 

R*8 

C2 

R*8 

DA 

R*8 

DE 

R*8 

DI 

R*8 

DP 

R*8 

DV 

R*B 

DX 

R*8 

DY 

R*8 

Dl 

R*8 

D2 

R*8 


Description 

semi -major axis 
variable used to compute v 
eccentricity 
inclination 
mean motion 
semi-latus rectum 
time 
V + w 

true anomaly 
X-coord, 

Y-coord. 

A ^ , intermediate parameter 
Ajs intermediate parameter 
Bjs intermediate parameter 
Bj » intermediate parameter 
C ] , intermediate parameter 
C 2 > intermediate parameter 
a, differential of semi -major axis 
e, differential of eccentricity 
i, differential of inclination 
pt differential of semi-latus rectum 
V, differential of true anomaly 
X, x-comp of velocity 
Ys y-comp of velocity 
Dj , intermediate parameter 
Dj , intermediate parameter 


42 


FORTRAN name Format 


E2 (200) R*8 

OM R*8 

WO R*8 

XO R*8 

DEL (200) R*8 

DOM R*8 

DWO R*8 

EPS R*8 

FJ2 R*8 

FMU R*8 

NEG 1*4 

NMN R*8 

NPN R*8 

NPO R*8 

NSS R*8 

OZT R*8 

RDA R*8 

RSQ R*8 

TXO R*8 

TYO R*8 

XOO R*8 

XOl R*8 

XIO R*8 

YOO R*8 

YOl R*8 

YIO R*8 

ALFO R*8 

ALFl R*8 

ALF2 R*8 

ALF3 R*8 

COST R*8 

COSU R*8 

COSV R*8 

COTI R*8 

DELM (200) R*8 

DOCI R*8 

DXOO R*8 

DXOl K*8 

DXIO R*8 

DYOO R*8 

DYOl R*8 


De scription 

dummy array of integrated values 
Q, longitude of ascending node 
Wq , argument of perigee 
Xq, distance between origins of 
dummy array of differentials 
^ differential of long, of ascending 
node 

Wqj differential of arg. of perigee 

€ 

J2O 

yu, - 3*986032 D5 

no. of parameters to integrate 

II - n 

n + n 

n + Wg 

n* sin2 i 

CO 

r^/a* 

r^, square of earth's radius 
intermediate parameter 
intermediate parameter 
Xqq, intermediate parameter 
X^j s intermediate parameter 
Xjg , intermediate parameter 
Y^g, intermediate parameter 
Yqji intermediate parameter 
Y » Intermediate parameter 
intermediate parameter 
a^t intermediate parameter 
a^i intermediate parameter 
intermediate parameter 
cos (i) intermediate parameter 
cos (u) Intermediate parameter 
cos (v) intermediate parameter 
cotan (i) intermediate parameter 
^ay of derivatives 
n . cos I, intermediate parameter 

?o, 

Too 


FORTTiAN Nam© 


Format 


Description 


DYIO 

R*8 

Yio 

ELEM (200) 

R*8 

array of integrated variables 

LA Ml 

R*8 


LA M2 

R*8 

LAMS 

R*8 


LAM4 

R*8 

\oi 

LAMS 

R*8 


LAM6 

R*8 

4° 

LAM7 

R*8 

LAMS 

R*8 

^4° 

LAMS 

R*8 


LAMIO 

R*8 

X 10 

k\o 

LAMll 

R+8 

LAM12 

R+8 

4° 

N.BAR 

R*8 

n 

SINI 

R*8 

sin (i) 

SINU 

R*8 

sin (u) 

SINV 

R*8 

sin (v) 

TANXJ 

R*8 

tan (u) 

TDXO 

R*8 

X (0), initial value of X-component 
of velocity 

TDYO 

R*8 

Y (0), initial value of Y -component 
of velocity 

TNMN 

R*8 

2n - n 

TNPN 

R*8 

2n + n 

ZERO (4) 

R*8 

array of initial values of position 
and velocity 

COS2I 

R*8 

cos (2i) 

COS2U 

R*8 

cos (2u) 

SIN2I 

R*8 

sin (2i) 

SIN2U 

R*8 

sin (2u) 
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ANALYTIC PROGRAM 


Genersd Overview 

The analytic program is comprised of a MAIN routine and four (4) sub- 
routines, The MAIN routine reads in the data (via 2 namelists) and writes out 
the initial conditions. Execution starts by integrating the equations of motion 
for the desired time jieriod. Values for the orbital elements are written for 
each time point. 


PliECEDXNG PAG£ 


blank not FILMi:^ 


47 




I 


4 




I 





X 



MAIN Houtine Description 

Purpose g To store initial values road in of the variables to be integrated. The 

variables not read in are calculated from those that are. The following 
are the variables that are integrated: 
a, 4, i, n* u, Wq , 
u is obtained from; 

U = V + Wq 

The integrator is then called. All values are calculated for each time point 
until the specified final time is reached. 
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Subroutine DERIV Description 

Purpose; To compute values of the differentials to be used by the integrator 
(FNOL2). 

Method; a, e, ijQj 14 Wq are calculated by equations (30) j (31) g respectively. 
The intermediate variables used to determine position and velocity 
are obtained from; 

Xq 

e 

n 


A |g 

B,,C, 

[ # 

A 2 * 

Bj.C 

2» 




^ 1 * 

A. 2* ^3 

9 * e « 

x„„ 


^ 0 - 

Too 

» Yol * 

T.oi 

Too 

» ?01» 

Tio 

Yoo 


Y,o 


k 


12 


These expressions are to be found in the listing of the analytical program in 
Appendix V, 
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Subroutine Description 


I, MAIN 

A, Purposei To read in control parameters and initial conditions. 
To write out control parameters and Initial conditions, 
Bo Csdls; FNOL2 

Co Method; All input through use of namelists; 


1, INIT; 


Variable 

X(0), X(0), Y(0), Y(0) 
a 
6 
i 
V 

n 



Variable Name 

Format 

ZERO (4) 

R*8 

A 

R*8 

E 

R*8 

I 

R*8 

V 

R*8 

OM 

R*8 

WO 

R*8 

FJ2 

R*8 


Description 

initial value of pos & vel, 
semi-major axis 
eccentricity 
inclination 
true anomaly 
longitude of node 
argument of perigee 
oblateness constant 


. INTEG 
Variable 

Variable Name 

Format 

Description 

At 

TF 

R*8 

final time 

TINCR 

R*8 

integration time step 


TI 

R*8 

initial time 

TNTP 

1+4 

print interval 


KUTMOL 

1*4 

integration type 


PRINT 

L*1 

flag for intermediate values 


NEQ 

1*4 

no, of variables to integrate 
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B. Common blocks useds 


1. ELEMENT 


Variable 

FORTRAN Name 

Format 

Description 

a 

A 

R*8 

semi'-major axis 

e 

E 

R*8 

eccentricity 

i 

I 

R*8 

inclination 

n 

OM 

R*8 

longitude of node 

V 

V 

R*8 

true anomaly 

u 

U 

R*8 


Wo 

WO 

R*8 

argument of perigee neglecting 
periodic behavior 


B. Common blocks used; 

lo ELEMENT 

A, E, I, OM, V, U, WO 

2. INTEG 

TF. TINCR, TI, INTP, KUTMOL 
3o CONST 

FJ2j FMU, RSQ, ZERO 

4. FLAG 

PRINT, J, NEQ 

5, A NS 

X, Y, XO, DX, DY 


n, DERIV 

A, Purposes To define variables to be integrated (differentials) 

B, Called by; FNOI2 

C, Method; Differentials defined in terms of variables previously defined 

D* COMMON blocks used; CONST; FJ2, FMU, RSQ, ZERO 

FLAG; PRINT, J, NEQ. 

ANSs X,Y,XO, DX, DY 
ALFA; ALFO, ALFI, ALF2, ALF3 
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APPENDIX IV 


c 


PURPOSE TO sntVE THE EOUATIUHS OF MOTION OF TWO SATELLITES NUMERICALLY 

FoF. COMPUTE VALUES flE POSITION! X,Y,/.»RI VELOCITY! K»»YSZSR“ >» 
SEMI-MAJOR AXIS! A )* ECCENTRICITY! £ ), INCLINATION! I I« 

LONGITUDE OF NODEt ARGUMENT OF PERIGEE» U,V AS A FUNCTION OE 
TIMF . 


PROGRAM WRITTEN FOR R. BARBIERI BY JOEL MASHRAUM 1-1-71 


INPUT 


RY NAMELIST 
1. NAMFLIST 
A. STEP = 
R. TI = 
C . T F = 
0, INTP = 
E» KUTflOL 
F. EPOCH 
G , J ?0 = 

NAMELIST 
A, X a 

R» Y 

C» Z = 

XOO? = 

YOOT = 


INTFG 

INTEGRATION STEP SIZE» TIME IN SEC, ( 600, I 

INITIAL time, secs ! 0.0 I 

FINAL TIME, SECS ! RA.AOO SEC = 1 DAY I 

PRINT frequency, integrations ! 1 * EACH STEP I 

= method of integration ! 2 = RUNGE-KUTTA fi AOAMS-MOULTON 

= GATE OF INITIAL CONDITIONS ! 720101.0 I 

OBLATENESS CONSTANT ! 0,0010823 I 

INI T 

X-CQORDINATE OF POSITION ! 1.0 I 
V-COORO»NATE OF POSITION ! 1,0 I 


Z-COORDINATE of POSITION ! 1,0 ) 
X-COOROINATE OE VELOCITY ! 0,1 I 
v-COORDIWATE of velocity ! 0,1 ) 


ZOOT 
G, A 
Ho E 


Z-COOROINATE OF VELOCITY ! 0,1 I 
SEMI-MAJOR AXIS ! 1500, KM ) 
ECCENTRICITY ! 0.01 ) 


I, I = INCL INAT!ON!DEG ) ! A5 . I 

J, OMEGA = LONGITUDE OF NODE, DEG I <>5, I 

K, W = ARGUMENT OF PERIGEE, DEG ! 0,0 I 
Lo V = TRUE ANOMALY, DEG ! 0,0 I 

H. omega's flag TO set WY £ WV® a 0 ! 1 I 

N, KPRINTs FLAG TO PRINT VARIOUS INTERMEDIATE VALUES ! 0 1 

O, KPLOT = FLAG TO DRAW 4060 GRAPH ! ,TRUE. I 


IMPLICIT REALMS !A-H,0-Z) 

REAL*fl L,M 

REAL*B r,J20 

LOGICAL KPLOT 

INTEGER PR, PI, YES, NO 

niMENr^ON ELEM!20)*DELM!20) 

dimension RI5000l,nR!5000) , T!5000I 

niMENSlON TYPEI5,3I 

COMMON /INCONP/ A , E , I *nMEGA ,W , U ,X ,V » Z , XDOT , YOOT , ZOOT 
COMMON /PLOTS/ T,R,OR,J 

COMMON /CONST/ TF , STEP , T 1 , I NTP , KUTMOL , EPOCH, FMU, J20»RSO 
COMMON /flag/ XO, KOMEGA,KPR!NT,KPLOT 
DATA ELFM.DELM / 40*0,00 / 

DATA YES, NO / "YES » , «N0 « /,PL,PR/ «N0 ','NO * / 

DATA TYPE /'RUNGE-KUS 'TTA THRO « , HiGHOUT »,2*’ S 

1 “RUNGE-KU* , "TTA AND S « ADAMS-MO • , ' ULTON » , ' * . 

2 'RIINGE-KUS'TTA WITH' ERROR C ‘ • ALCULAT 1 » , " ON */ 

NAMFLIST /INTEG/ STEP , T I , TF , 1 NTP, KUTMOL, EPOCH, J20 

NAMELIST /IN IT/ X,Y,Z ,XOOT ,Y [XiT , ZDOT, A, F, I , OMEGA , W, V.KOMEGA , KPR INT 
, , KPLOT 

RAOIAW(XI = X / 57,2957795130823 
REA0!5,INTFGI 


REA0!5,INITJ 
I a RAOI AN! I I 
W = RAOIANtW) 


PBECEDING PAGE BLANK NOT FIUl 
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V » RADIA(^«Vt 

OMEGA s RADIAN(OMEGA) 

L ■ E » OCOS(W[ 

M ■ E • nSIN(Wf 

U»W+V 00005700 

ELEM( U « X 

ELEMl 2) a V 

ELEMJ 3» a Z 

FLEMJ <») a A 

ELEMJ 5( a L 

ELEMJ 6) a I 

ELEMJ 7t a OMEGA 

ELEMJ 8) a M 

ELEMJ 9) a U 

ELEMJ lOJaXDOT 

ELENJllJaVOOT 

ELEMJ12PaZ00T 

LL * 0 

WE a 0 

WRITEJ6«U 

1 FORMAT { *1“ ) 

WR|TEJ6*2J 

2 FORMATJ 30X» lOJ 1H*( ♦ 10K» ‘ INITE AL COWOITl nWS' V IOX»10J IH*) » 

CALL OATEJ EPOCH J 

WAITEJ6»3> X*Y»Z 

3 FORMATJ •-POSITION VECTOR SlOXc'X a • »01 5.5. 10X» • V = S015.5. 

U 10X.«Z a ’.OIS.S » 

WRITEJ6.AI XDOT, YOOT, ZOOT 

A FORMATJ •OVELOCITY VECTOR »,9X.«XM = • , D1 5 . 5 .9X . • V • • » *.015.5. 

U 10X.«Z»« a »,015.5 ) 

WRITE J 6. 5 I A.E.I .OHEGA.W.V.JZO.KOMEGA 

5 FORMATJ /39X, “SEMI-MAJOR AXIS JAI * • .020. 10/A2X . • ECCEWTR IC I TV JE» 

1 a • .020. 10/A2X." INCLINATION JFII * *. 020, 10/23X ,* LONGITUDE OF*. 

2 • ASCENDING NODE IOMEGA) a • ,020. 10/35X, • ARGUMENT OF PERIGEE JW) 

3« • .020, 10/A2X. 'TRUE ANOMALY JV) = • .D20,10/«2X. “ORLATENESS JJ20) 

4a • »n20,10/52X. “OMEOAX a 8,I2 ) 

IFJ KPLOT ) PL a YES 

IFJ KPRINT.NE.O ) PR a YES 

WRITEJ&.6) PL.PR 

6 FORMATJ ////45X. “PLOT = • . A4.5X . «PRINT a • .A4 ) 

WRITE J6, 7) JTYPEIJ.KUTM0L)*Ja|.5).STEP.INTP 

7 FORMAT! ////48X . « I NTEGRAT ION PAR AMETERS • /48X » 22J IH- ) //45X. 

1 “METHOf) SSAS/ASX.'STEPSIZE • . 020. 10M5X , • PR INT FREOUENCV 

2 .IR ) 

WR f TE f 6 9 i ) 

call FN0L2J KUTMOL.IZ.STEP.LL.INTP.NE.TI.ELEM.OELM) 

IFJ KPLOT ) CALL PLOT 
WR1TEJ6.9) 

9 FORMATJ »-BNO OF RUNM 
RETURN 
ENO 




SUBROUTINE HATE ( EPOCH » 

IMPLICIT REAL*8 »A-H»0-Z) 

INTEGERPA YR*DAeHReSEC®VEAR 
DIMENSION MONTHUai 

DATA MONTH/ • JANN » FEBS ® MAR ' * » APRS » MAYS “JUNES 
U “JULYS > AUGS “SEPTS * OCTS ' NOV" v “ DEC“ / 

YR » EPOCH / SD^> 

YEAR « 1900 + YR 

Jl ■ EPOCH - YR * 1,04^ 

MO « J1 / 1.02 
OA » Jl - MO • 1,02 

HI ■ EPOCH = YR ♦ l.D^ - MO * 1.D2 - DA 

MR » HI • 1.02 

MIN ■ HI * l.DA ~ HR * 1.02 

SEC » HI • 1.06 - MIN * 1.D2 - HR • I, DA 

WRITE 16,51 MONTH! MO S DA, YEAR, HR, MIN, SEC 

5 FORMAT! // //AAX , “EPOCH » “ , AA, I 3, “ , “ , I A ,5X , 1 2, “ HRS “,I2,“ MINS 
$ ,I2*“ SECS“ ) 

RETURN 

END 
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SKRRriHTINf? 'ERrv(T,ELEM,OELMI 
implicit RF AL*fl( A-H,n- 2 I 
REAL*fi L«M 
Lnr.iCAL *SFN»icns 


VERSION OF 
SUBROUTINE 

INPUT 

T 

FLFM 

OUTPUT 

OFLM 


12-23-70 

DERIV COMPUTES THE DERIVATIVES OF THE ELEMENTS ARRAY, 


-TIME 
- ARRAY OF 


ELEMENTS 


ARRAY OF DERIVATIVES OF ELEMENTS 


00000 1 00 
00000200 

00000300 

oooooAon 

00000500 

00000600 

00000700 

OOOOOROO 

00000900 

00001000 

00001100 

00001200 


rLFM 

DELM 


0000 1 300 

X 


1 

00001 AOO 

Y 

V « 

2 

00001500 

Z 

1 « 

3 

00001600 

AO 

A0» 

4 

OOOOlflOO 

t.O 

LO' 

5 


10 

I0«» 

6 

00002000 

OMf:OA 

OMEGA* 

7 

00002100 

MO 

Mf)< 

8 


UO 

U0» 

9 


x« 

X s * 

10 

00002A00 

V • 

Y 8 fl 

11 

00002500 

1 « 

Z » ® 

12 

00002600 


DIMENSION FLEMI 11 ,OELMn ) ,ELF( 20) ,0EL(?O> 

COMMON /FLAG/ XO, KOMEGA »KPRINT , KPLOT 

COMMON /CONST/ TF » STE P , T I , INTP »KUTMOL, EPOCH , FMU, F J7 , R SO 
EQUIVALENCE lELE n)»X),(6LE (?),Y),(ELE (3)»Z), 

UELE (A)»A),(ELF (5),L)9(ELE (6l,FI),(ELE (7),DMI,(ELE (8)»M)» 

21 ELF « 9) »U) » IDEL m »DX,ELE( 10) 1 ♦ (DELI2) »0Y ),(0EL(31 ,D2 ) * 

3IDEL (A)»OA),(DEL (Bl^DLlsIDEL (6)»DI)p(DeL ( 7 ) » DOM ) , ( DEL (S),DM) 
A (DEL I 9) 9 DU) ,(DEL< 10) gDOX ) , (DELI 11 ) *DDV ), (061(12) »0D2) 

RESTORE INPUT ARRAYS FOR EQUIVALENCE STATEMENT 
00 20 I*1«10 
ELFd )=6LEM( I ) 

20 OEU n=DELM( I ) 

NAMELIST /TRIG/ COSV, SINV»COSU»SINU,COS! » SI NI * C0S2U* S IN2U» S 1 N2 I , 

U C0S2I ♦TANU»COTI 


A2»A3»B2»B39C2»C3 
V»UtF I >P«XO*FN 

UA»DL»DP»C» 0X0»DF ,DOM,DM,OUtDOI 
WX fMV «WZ ;DWY fOWZ 
nx,0Y,02»DUX»D0Y9D0Z 


20 


2001 

2000 

500 


NAMELIST /NAMl/ 

NAMELIST /NAM2/ 

NAMELIST /NAM3/ 

NAMELIST /NAMAZ 
NAMELIST /NAM5/ 

DO 20 ! = U12 
ELFd )»6LEM( I ) 

IFF KPRINToEO.O 
WRITE(6,2001 F 
FORMAN so ENTERS ) 

WRI TE( 6 » 20 nO) )ELEM( I ),0ELM( I )d*l»13) 
22X* sflEMENTS' 99X»sdERI VATIVESS //( 10X,2020o8) ) 


00002700 

00002ROO 


00003500 

00003600 

»00003fl00 

OOOOA200 

OOOOAAOO 

0000A500 

0000A600 

OOOOA7UO 

0OOOA8OO 

OOOOA900 

00005100 
00005200 
00005300 
0000 5 AOO 


) GOTO 500 


FORMAT! 

CONTINUE 
6 » DSORT( L*L + M*M ) 
IFF DARS(M) ,LF. 1,0-30 
W = 0 ATAN 2 F M«L ) 

GOTO 60 

io w a n,no 

60 V » U - W 


oANO. DAflSF L) ,LEol ,D-30 ) 


GOTO 50 


COMPUTE FREFJUENTLV USED TRIGONOMETRIC FUNCTIONS 


00005600 
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cosv«ocns(v ) 

SINV>nSIN 9 V) 

C 0 SI 1 «DC 0 SIU» 

C 052 U>DC 0 S 12 .D 0 *(»» 

SIN(i«DSINIU) 

SIN 2 UBn 5 IN( 2 . 00 * 0 ) 

SINI>OSIN(PI) 

$t^( 2 l« 0 SIN« 2 oD 04 'Fn 

cosi>ncosiFt » 

COS 2 r*DCOS« 2 .DO*FI ) 

TANU>OTAN(Ut 
COTI«OCnTAN*F| ) 

IFtKPRINT.eO.U HRITEt 6 «TRIGt 

P«A*a . 00 “F*FJ 

FN«DS 0 RT(FMU/A** 3 J 

DENOM * 1.00 + L * COSH + M * SINU 

XO « A l.DO - L*L - t / DENOM 

IFlKPRINT.EO.n HRnE< 6 sNAM 2 ( 

OA«- 3 .nO*FMU*Fj 2 *RSQ/«FP<*OSORm.DO-E*En*( ( L* SINO - M • COSU ) 
• ♦{ 1,00 - 3,00 • S!NI ♦ 

1 S!NI*S!NU*SINUI/X 0 ** 4 +P *SINI*SIN 1 *SIN 2 U/X 0 ** 5 ) 

DI«-n , 500 *FMU*F J 2 *RSO*SIN 2 I*SINU*COSUf /«FN*A*A*DSORT« 1 ,D 0 -E*EI 
1 ♦X 0 ** 3 ( 

OOM« 01 *TANU/«SINI» 

OL ■ M • DOM * COSI - 1,500 ♦ FMU ♦ FJ 2 * DSORTl 1.00 - E*E i* RSO 

1 / FN / A SINO •( 1,00 - 3,00 ♦ SINI ♦ SINI>t<SlNU * SINU 1 / 

2 XO*#A +( XO* L / P +« 1 DO XO / P (* COSU 1 * SINI * SINI 

3 SIM 2 U / X 0**4 1 

DM = -L ♦ DOM * COSi + 1.500 ♦ FMU ♦ FJ 2 ♦ RSO ♦DSORTt 1.00 - E*E 

1 / FM / A COSU l.DO - 3,00 * SiNl ♦ SIN! * SINU ♦ SINU )/ 

2 X 0 **<» -( XO* M / P +1 l.DO + XO / P f* SINU 1 * SINI ♦ SINI ♦ 

3 SIN 2 U / X 0**4 I 

DP - DA *( l.DO - L*L - M*M » - 2,00 * A *( L<-OL + M*DM ) 

OU » FN • A * A * DSORTl l.DO - E*E I / X 0**2 - DOM • COSI 

0 X 0 = ( OA <■( 1.00 - L*L - K*M 1 - 2,00 * A L ♦ DL + M • DM M 
t DENOM - A •( 1,00 - L*L - M*M )*( OL ♦ COSU - L • DU * SINU 

C + DM • SINU + M ♦ OU • COSU 1 / DENOM / OENOM 

IFIKPRINT.EO.n WRITE( 6 »NAM 3 t 
W Z«DSORT ( FMU*P 1 /X 0 *A 2 
IFIKOMEGA.EO.OI GO TO 30 
WX«D|/COSU 
WVsO.DO 

DOI = (-l . 5 OO*FJ 2 *RS 0 * 0 S 0 RT(FMU/P»/XO** 3 >* U DA-A* 0 P/ { 2 , DO*P » 

1 - 3 . 00 *A*DX 0 /X 01 *S IN 2 I*SI NU*COSU+A*(D! * 00521 *S 1 N 2 U 

2 +OU*SIN 2 I*COS 2 un 

DW X» I noi * 0050 + 01 *OU*SINU 1 /COSU*COSU 
DWY« 0 .D 0 
GO TO <£^0 
30 WXrO.DO 
WY-O.DO 
OWX-O.DO 
OHY* 0.00 

AO DMZ*M?*(OP/ 12 . 0 O*Pl- 2 , 00 *OXO/X 0 ) 

IFIKPRINT.EO.U WR I TE t 6 ,NAMA I 

A 2 a 2 , 00 *(MY*DZ-DV*W 2 J+ 2 * 0 MV-Y*DMZ+WX*(Z*WZ + Y*WY» 

A 3 =-WZ*WZ-WV*WV + FMU/X 0 ** 3 -( l, 5 D 0 *FMU*FJ 2 *R S 0 /X 0 ** 5 f 
1 *( 1 ,D 0 - 3 ,D 0 *SINI*SINI*SINU*SINU) 

B 2 « 2 .D 0 *« 0 X*WZ- 0 Z*WX!+X*DWZ-Z*DWX+WY*IZ*WZ+WX*! X+XOl 1 
B 3 »-«X*«X-MZ*MZ+FMU/X 0 ** 3 -a, 50 O*FMU*FJ 2 *RS 0 /X 0 ** 51 *SINI*S!NI* 
i C 0 S 2 U 

C 2 « 2 , 00 *IOY*WX-MV*SDX+DXO! 1 +V* 0 MX-X*DWY+WZ* «WY*V+WX*X)-XO*OWY 

C 3 --MV*«V-NX*WX+FMU/X 0 ** 3 -( l. 500 *FMU*FJ 2 *R S 0 /X 0 ** 5 !*SINU*C 0 S 2 I 

IFIKPRINT.EO.n WRITEI 6 ,NAMn 

D 0 X«-A 2 -A 3 *X 

0 DY»-B 2 -B 3 *V 

OOZ«-C 2 -C 3 * 2 -XO*WX*WZ 


00005800 

00005900 

00006000 

00006100 

0CCC6200 

00006300 

00006A00 

00006500 

00006600 

00006700 

00006800 

00006900 

00007100 

00007300 


00007600 

00008600 


00008900 

00009500 


00009100 

00009200 

00009300 


000 10 100 
00010200 
00010300 
OOOIOAOO 
00010500 
00010700 

00010800 

000 1 1000 
00011100 
0001120 

00011300 


00011600 

00011700 

00011800 
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OOOH90O 

00012100 


100 DO 110 I»U12 
110 oeiMt n«DEi.n > 

IP(KPRI»«T.E0.n WRITe(6«NAM5) 

IFC KPRINT.eOoO » GOTO 99 
MRITE(6»2002> 

2002 FORMAT MO EXIT’ ) 

WRtTE<«S,2000t 4ELFM( n«0ELM(I)9l>ia3) 

99 RETURN 
END 


00012500 
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SUBROUTINE PLOT 
REftL*8 T!M6»POS»VEL 

DIMENSrON TIMEtSOO) » VEL ( 500 > t POS ( 500 > 
dimension T(500)»RI500»9V<500».Z(200) 

COMMON /PLOTS/ TIME*POS»VEL»J 

C CONVERT FROM DOUBLE PRECISION ARRAYS TO SINGLE PRECISION 
00 10 1= 1 »J 
Tin * TIME ( n 
R«n * posi n 
10 V« n = VEL(I ) 

CALL MODESGI Z«0 ) 

CALL SETSMGI Z»84»'+“ ) 

CALL GRAPHGt Z 9 J »T «R , 13» » T 1 ME ( SECS M 9l5*"POSITION ( KM M917, 

U 'TIME VS« POSITION* I 
CALL LINESG( ZjJoT.R ) 

CALL PAGEG( Z»09l«l ) 

CALL GRAPHS! Z « J ,T 9 V» 13» * TI ME ( SECS I « , 19 , « VELOC IT V ( KM/SEC IS 
U 17, 'TIME VSo VELOCITY' ) 

CALL LINESGt Z 9 J 9 T 9 V I 
CALL PAGEG( Z,0«1»I ) 

CALL EKITG ! Z I 

RETURN 

ENO 
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SUBROUTINE EN01.2 ( J«N«G * L *M»NE « X «Y « D ) 

IMPLICIT REAL*B( A-H,0-Z,$) 4480 

4500 

4510 

DOUBLE PRECISION X0,V0,VA,VC,VP,Y1 4520 

DIMENSION Y(50)»n(50>,YB(30o6)«GI2I30l«r,I 3I30I«GI4(30) «EF( 30) , 4530 

IEF1(30),EF2(30) .EP3(3O)tYl(3O)«ERRORI30l eHA(30),YAI50l .0A(50I« 4560 


2VC(30)*VPI30),Y0I50I 
EC»Y«N+3) 

1 HaG 

2 HZ»H 

3 LN»N-«-MAX0(L»3> 

SUBROUTINE FN0L2( JfN*G«L«M,NEvX,Y 

4 NA»0 

5 NB«1 

6 NF»0 

7 NG=0 

8 FaOoOO 

9 FAaO,00 

10 FB*0o00 

11 Fc«o.no 

12 CONTINUE 

13 ENEaNE 
00 200 1=1 ,LN 

200 VDm*DBLE(V( I I I 

xd-obleixi 

DO 200 I=1»LN 

200 VOI n-YU) 

XD=X 

14 IFIJ-3I15, 21*15 

15 IF(NE)18*16*18 

16 JA=4 

17 GO TO 22 

18 Rei«10.00*4<-eN6 ) 

19 RE2=10.004«(«ENE-3.000) 

20 R6M«10.00**I-ENE-1.500) 

21 JA=1 

22 00 25 I=1*N 

23 DO 24 IC*1*5 

24 VBn*IC)«0.00 

25 ERROR (I >=0.00 

26 CALL DERIV(X>Y*D) 

CALL TERMtXtY*0»F) 
!F(DABS(F)-l.00o9) 731*731*5209 

5209 CONTINUE 

DO 300 1*1 »N 
GI2U )*om 
Gi3(n*Dm 
GI4(I >*D(1I 
300 EF(n*Dm 

C 27 CALL OUTPUTfX»Y*D*ERROR*N*L*H) 

27 CALL OUTr:»V*D*ERROR*N*L*M» 

28 FD=Y(N+1) 

29 IFIJ-2I 30*129*30 

30 GO 70(31*37*35*37), JA 

31 DO 33 l»l*LN 

32 YA( I )*V0( I ) 

33 DAmaOin 

34 GO TO 37 

35 HB*H 

36 M*2.00*H 

37 H02 « .500*H 
00 39 |«1,N 

38 YB(I*N8l«0m 
XL • 0(1) « H02 

C 39 Ym«SNGL(YD(I)«XL) 


4570 

4580 

4600 

4610 

4620 

DERI V, TERM, OUT) 4630 

4650 

4660 

4670 

4680 

4690 

4700 

4710 

4720 


4750 

4760 

4770 

4780 

4790 

4800 

4810 

4820 

4840 

4830 

4850 

4860 

4870 

4880 

4890 

4900 

4910 

4920 

4930 

4940 

4950 

4960 

4970 

4900 

4990 

5000 

5010 

5020 

5030 

5050 

5060 

5070 

5080 

5090 

5100 

5110 

5120 

5130 

5140 

5150 

5160 

5170 

5180 
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39 Yn»»VD(H*XL 
C X«SNGL(XD+H02» 

X>XD4-HD? 

40 CALL OERIV lX»V9GI2t 

41 DO 42 I«1«N 

XL * GI2(n*H02 
C 42 YU >=SNGLiVO<l )+XL( 

42 V<n-Y0nKXL 

43 CALL DERI9 <X«V«G131 

44 00 45 I»1,N 
XL»GI3n )*H 

C 45 Vm«$N6L(Y0n ( + XL> 

45 Y(n«YOm+XL 
C XaSNGLfXO^H) 

X»XO+H 

46 CALL 0ERIV(XyY«GI4» 

47 H06 »H/6o00 

GO Tat48»55*60»66> »JA 

48 DO 52 l«lvN 

XL«(D(IP ♦ 2.00*(GI2n> + GI3n»» +G14UJ>*H06 

49 VC< I )*YD( n+XL 

51 YD<U«YAnJ 

52 ERROR I n *0.00 

53 JA*3 

54 GO TO 35 

55 DO 57 !»1,N 

XL»(D(n ♦ 2.D0*»G!2(H + GI3tn» +GI4(n>*HD6 

56 V0< n=VD(n+XL 

C 57 ERR0Rn>»SMGUY0m-VPnP)/15, 

57 ERROR( n»(YD(n-VPn ) P/lSoDO 

58 JA»1 

59 GO TO 681 

60 00 62 1*1 

61 vo( n=Ycn ) 

XL=i(Om + 2,00*(G12<n + GI3(m 4GI4m>«HD6 

62 VPm»VAnJ+XL 

63 K«HB 

64 JA*2 

65 GO TO 681 

66 DO 68 I«t*N 

xL«inm + 2o00*(6i2m cnnu •^6i4(n)*H06 

67 vom»vom+xL 

66 ERROR (I MO. DO 
661 DO 69 I*1«N 
C 69 vn MSNGLtYOn ) ) 

69 vm«Yom 

XD*XD-i-H 

C X*S»«GLtXD) 

X*XD 

70 CALL DERIV«X.Y»OP 

71 FC*F 

72 CALL TERM«X,V^O»Fl 

73 IF(0ABS«F)-l.0D-9 (73197319733 

731 NF*5 

732 GO TO 124 

733 IF1FP74»124976 

74 FA*1.00 

75 GO TO 77 

76 FB«1.D0 

77 IF(FA"FB)S3fT8t63 

78 MF*NP'*’l 

79 JA>4 

80 NB»1 

81 H«H*F/«FC-FJ 

82 IP{NF<>4137»37,I24 

83 !FINE184«117984 
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5190 

5200 

5210 

5220 

5230 

5240 

5250 

5260 

527 

5280 

5290 

5300 

5310 

5320 

5330 

534 

5350 

5360 

5370 

5380 

5390 

5400 

5410 

5420 

5430 

5440 

5450 

5460 

5470 

5480 

5490 

5500 

5510 

5520 

5530 

5540 

5550 

5560 

5570 

5580 

5590 

5600 

5610 

5620 

5630 

5640 

5650 

5660 

5670 

5680 

5690 

5700 

5710 

5720 

5730 

5740 

5750 

5760 

5770 

5780 

5790 

5800 

5810 



5630 

5840 



84 IF( JA-1 in7«8S» 117 S8S0 

85 IFIJ“3)86,117*86 5860 

86 DO 95 I«l*N 5870 

IFIVnnS86«885«886 5880 

885 HAU) >1000^00 5890 

GO TO 95 5900 

886 IP(EC »880«390«87 5910 

87 IFIDABSIVt I n°eC > 880,880*890 5920 

880 IF(0AB$(ERR0R(1 MH1E2) 882,94,881 5930 

881 IF(DA8S(ERR0R(I n>REl>94,94,882 5940 

882 HAm REM/ (DABS! ERROR n H + . 0000000001 DON **(.200) 5950 

863 GO TO 95 5960 

890 IF(nAeS(ERRORm/Vn n-RE2)892,94,B9l 5970 

891 tPIDABSlERRORn )/vn n-REl)94,94,892 5980 

892 HA( I )«H*(REM/(OA8S(ERRORn )/T(I 000000000100))**! .200) 5990 

893 GO TO 95 6000 

94 HA( I )«H 6010 

95 CONTINUE 6020 

96 HB«OA0S(HA«NH 6030 

97 DO 98 I«1»N 6040 

C 98 H6»AMlNHABS9HAn ) ),HB) 6050 

98 HB>DMlNnOAB5(HAn n ,HB) 6060 

99 !F(OABS(H)-MB)100,117,10l 6070 

100 IF(0ABS4HZ)-0ABS(HH101,101,116 6080 

lot DO 103 l-l,LN 6090 

102 YD(n>VAm 6100 

C VI I)»5NGL(V0tl n 6110 

Y(I)«YDm 6120 

103 om-DAd) 6130 

104 IF(NB-6) 107,105,105 6140 

105 XO-XO-H 6150 

106 60 TO 109 6160 

107 X0«X0“2.D0*H 6170 

108 HZ > H 

109 H*DSIGN(HB,H) 6190 

X«XD 6200 

C X>SNGL(XO) 6210 

CALL 0ER!V(X,Y,D> 6220 

CALL TERM(X,Y,0,F) 6230 

no NB»1 6240 

111 XABS«OABS(.000001DO*X) 6250 

112 IF(DA8S(H)«XA6S)113,n3,117 6260 

113 NG«NG*1 6270 

114 H»DSIGNiXA6S,H) 6280 

115 IF(NG > 10)124,150,150 6290 

150 MRITE(6,1261) H 6300 

1261 FORMAT ( IHI ,107HEXECUT]nN TERMINATED BECAUSE INTERVAL OF INTEGRATI 6310 
ION LESS THAN l.OE "6 TIMES INDEPENDENT VARIABLE (X). X a,lPD15.7) 6320 

STOP 6330 

116 HZ*H 6340 

117 IF(M)n8,ll8,121 6350 

118 IFI (Y8N+1 )-FD)-V(N+2) )29,119,ll9 6360 

119 F0*FD+V(N+2) 6370 

120 GO TO 124 6380 

121 NA-NA*1 6390 

122 IF(M>NA)123,123,29 6400 

123 NA»0 6410 

’.124 CALL OUTPUT(X,Y,D, ERROR, N,L*H) 6420 

124 CALL OUT(X,Y,0, ERROR, N«L«H) 6430 

125 IFINF"4I29*29,126 6440 

126 HRITE (6,127) 6450 

127 FORMAT (IHO) 6460 

128 RETURN 6470 

129 NB«NB*| 6480 

130 lF(NBo6l30,l31,136 6490 

131 DO 134 tai,N 6500 
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132 EF3H J-VBUf,3l 6510 

133 EF2m-V0n*4> 6520 

134 EFUI )>Yfin«5S 6530 

135 GO TO 137 6540 

136 NB*10 6550 

137 H024 *H/24.D0 a560 

00 13S l>l«N 6570 

XL »(55.00*Dm -59,D0*EFim 437,00*EF2n ) -9,004EF3< IN *H024 o:>«0 

YPm«VDm4XL 6590 

138 Yl n>$MGUYPII n 6600 

X>SNGLIXD4-H) 6610 

138 Y(n«YPm 6620 

X>X04-H 6630 

139 CALL OEPIVUtY*EFI 66 

140 00 142 !"1«LN 6650 

141 YA< S 1>YDI 1 1 6660 

142 DA( n-0( n 6670 

143 DO 148 I«l9N 6680 

XL • (9,00*EFllt ♦I9,00*0m -5oOO*EFUn '^EF2 ( | H *HD24 6690 

144 V0( n»VDU )4XL 6700 

C 145 ERROR n l»-SNGL(VDm-YP( n [/14. 6710 

145 ERRORd !»“(YOn J-VPU n/14.D0 6720 

146 EF3U »-EF2m 6730 

147 EF2(n-EFint 6740 

148 EFl n »«0( n 6750 

999 CCriTINUE 6760 

149 GO TO 681 6770 

END 6780 


SUBROUTINE TERNt T « E LEN , LM» COND ) 

IHPLICIT RFAL*0U-H,Q-2 ,$ J 

COMMON /CONST/ TF ,STEP*TI » INTP,KUTMOLeEPOCHsFMU*FJ?»RSO 

DIMENSION ELEMUK OELMtU 

CONO r TF - T 

RETURN 

END 


I 


% 


J, 




SUBaoUTtNE OUT tX«V 
IMPLICIT REAL* 84 A-H, 0 - 2 » 

DIMENSION T< 9 »» 049 ),£RR 0 R( 2 ) 

DIMEMSIOM RfSOOOIt V 45000 U T( 5000 > 

DIMENSION OEG 4 S) 

COMMON /FLAG/ XO* KOMEGA 9 KPR I NT* KPLOT 
COMMON /PLOTS/ T»R*V»J 
DATA TPI/ 6.2831853072 / 

VMAGU«Y»Z) « DSORTt X«X * V*Y * 1*2 I 
DE6REE4XI * X * 57.2957795130023 
J ■ J + I 

R 4 J> > VMAG 4 V( UvY( 2 )«Y( 3 > » 

VIJI«VMAG( D(UfDt 2 )«D( 3 t ) 

T«JI ■ X 
WMTE 46 « 3 I 
3 F 0 RMAT 4130 UH-M 

HRITE 46 « 4 I RUI.V(J) 

^ FORMAT! » 0 S 25 X, 'POSITION - • , 020 . 10 . 15 X, 'VELOCITY » ' ,D 20 . 10 } 

E * VMAGI Y 45 I.Y( 8 I«O.DO I 

IF! DA 85 (Y« 5 H.LE. 1 . 0-30 .AND. DABStY 4 8 H .LE. 1 . 0-30 I GOTO 50 
H B OATAN 2 4 Y 4 SI.Y 45 I 1 

GOTO 60 
50 W ■ 0.00 
60 CONTINUE 

V2 * Y49I ~ W 

V 2 ■ OMOD 4 V 2 .TPI I 

W > DM 004 H.TP! ) 

Y 46 I • DM 004 Y 46 I.TPI I 

V 47 J • 0 M 0 D 4 Y 47 I.TP 1 J 

V 491 * DM 004 Y 49 >,TPI » 

OEGm B DEGREE! Y 4 7 » I 

OEG 421 > DEGREE! M I 

DEG ! 31 « DEGREE! Y! 9) ) 

0 EG 4 AI B DEGREE! Y( 6 H 

DEGI 5 I B DEGREE! V 2 ) 

WR!TE! 6 . 100 n X . ! Y 4 II . I b l, 39 

1001 FORMAT! '-'.AX, 'TIME (SEC) » ' , 015 . 8 , 8 X. 'X 4 KMI * ', 01 S.S, 8 X, 

1 'Y 4 KMI B • .D 15 . 8 , 8 X .'2 !KMI b 015,8 I 
WRITE 46 . 1002 ) XO. ( 04 1 ) , I » 1 , 3 ) 

1002 FORMAT! 'O', 7 X,'X 0 !KM) b ' , D 15 . 8 , 3 X , • X " (KM/SEC) « '.DIS.B.SX, 

2 «Y" 4 KM/SEC) • • ,D 15 , 8 . 3 X,'Z" 4 KM/SEC) = ', 015.81 
WRITE ( 6,10031 Y 4 Al.e.V( 6 l. 0 EG(A) 

1003 FORMAT! ' 0 ' , OX, 'A IKMI - ' ,D 15 . 8 , 9 X , 'ECCEN > ' , 015 . 8 , lOX INCL « 

3 ,D 15 . 8 ,' RAO « ', 015 . 8 ,' OEG'I 
WRITE! 6 , 100 A) Y 47 ),OEG!l), W , 0 E 6 ! 2 ) 

1004 FORMAT! ' 0 ', 8 X, 'L.A.N. » ', 015 . 8 ,' RAO * ', 015 . 8 ,' OEG'.lOX, 

4 'ARG OF PER » ', 015 . 8 ,' RAD « ', 015 . 8 ,' OEG'I 
WRITE 46 , 1005 ) V 2 ,DEG (3 I ,Y 49 I , 0 EG 4 31 

1005 FORMAT! ' 0 ' , 13 X, «V b %D 15 . 8 ,' RAO b ', 015 . 8 ,' 0 EG', 19 X,'U b 

5 015 . 8 ,' RAO B ', 015 . 8 ,' OEG'I 
WR|TE 46 *i 006 l Y 45 I,Y 46 I 

1006 FORMAT! 'O', 13 X, 'L » ' 015 . 8 , 45 X , 'M ■ ' 015.81 
RETURN 

END 
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BLOCK DATA 

IMPLICIT REAL* 8 U<»MtO-Z 
REAL «8 I 
LOGICAL KPLOT 

DIMENSION R(5000)t VI5000)« T(SOOO) 

COMMON /ELAG/ XO» KOMEGAfKPR I NT«KPLOT 
COMMON /INCONO/ A«E « I .OMEGA»W «U vX «Y « 2 « XDOT« VOOT» 200T 
COMMON /CONST/ TP » STEP » Tl « INTP «KUTMOL«EPOCH«f MU*P J2«RSO 
COMMON /PLOTS/ T*R»V«J 

DATA PMU 9 RSO / 3*986032D5« 4*06809807/ 

DATA T*R,V*J / 18000*0*00*0 / 

DATA A*E*!»nMEGA*0 *M/ 1 .504* i .D->2« 3«45. DO* 0.00 / 

DATA STEP«TI»TP*INTP*KUTMOL/ 6.0?.0.00*R,6404* 1 *2 / 

DATA EPOCH] 720101.00 / 

DATA X»*« ,2 / 3*1.00 / 

DATA XDOT*VOOT*ZOOT / 3*1.0-!/ 

DATA PJ2 / 1.08230-3 / 

DATA KOMEGA.KPRINT, KPLOT / 1*0*. TRUE. / 

END 
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APPENDIX V 


TWO SATELLITE MOTION STUDY - ANALYTIC THEORY 
LAST REVISION 10-16-70 

implicit realms IA-H,0-Z) 

REAL*S I 
LOGICAL PRINT 
DIMENSION TYPE<5«3> 

DIMENSION 2FR0(4) «ELEM(20I fDELM(70I 
COMMON t ELMENT f A »E » I «OM«V ,U«WO 
COMMON /INTEG/ TF 9 T INCR »TI , 1 NTP« KUTMOL 
COMMON / CONST / F J2,FMU9RSO«ZERO 
COMMON /FLAG/ PRINT tJ«NEO 
COMMON /ANS/ X«Y «XO»DX«OY 
RAOIANIXI « X / 57.2957795130023 

DATA TYPE /"RUNGE-KU S " TTA THRO « , » OGHOU T *f2** •* 

1 •RUNGe-KU",»TTA AND * ♦ • AOAMS-MO® , » ULTON • , • S 

2 *RllNGE-e<U>»«TTA MITH«»» ERROR C * * ' ALCULAT 1 S *0N ' / 

NAMELIST /INIT/ ZERO * A« E » I 9 V vOM* W0« F J2 

NAMELIST /INTEG/ TF « T INCR « TI » INTPfKUTMOL* PRINT 9 NEO 
REA0(5«1NIT> 

REAOIS.INTEGt 
WRITE ( 6 ,U 
HRITE(69l001) 

1001 FORMAT! 25< / ) »35X . • ANALYTI C SOLUTION OF TWO-SATELLITE NOTION * 

F 9« STUDY ® 1 
WRIT 6 < 6 *n 
I FORMAT! »l« J 
WRITEfGfS) 

5 FORMAT! 30X.10! 1H« > t lOX « ® IN I TI AL CONDI T IONS® » IOX 9 10! IH* ) ) 

WRITEI6«0I A9E«l90M9W0 »V 9 FU 29 ZEROI 1 1 9 Z ERO ! 3 I 9 ZERO! 2 > 9 ZEROIA ) 

8 FORMAT! / 39 X 9 *SEMI-MAJOR AXIS !A) = S020, IO/A 2 X 9 * ECCENTR IC I TY !E> 

1 » ' 9020 . 10 / 42 X 9 ‘INCLINATION ! I) « ® 9D20ol0/26X* 'LONGITUDE OF', 

2 ' ASCENDING NOOF ! 0 M) s D20. 10/ 34X , ® ARGUMENT OF PERIGEE !W0I «® 

3 021,10/ 42X9'TRUE ANOMALY !VI » * 9020.10/42X, ®0BLATENESS U20I 
4« • , 020 , 10 / 57 X 9 'X =®02'3,10/57X, ®Y »' 020 . 10/56X, ' X " «'020.10/ 

5 56X,'Y" =',020.10 1 

WRITE 16 ,71 !TYPE! J9KUTM0LI 9 J=l95I,TINrR 9 INTP 9 NEO 
7 FORMAT! ////48X , ® INTEGRATION PARAMETERS® /48X922! IM-I //45X, 

1 ®METMOn ®95Afi/45X,®STEPSIZE ® ,020. 10/45X9'PRINT FREQUENCY 

2 9I8/45X, 'NUMBER OF EQUATIONS = ®I2» 

WRITF!6,U 

U = V ♦ W0 
I « RADIAN!! I 
OM « RADIANIOM) 

U « RAO I AN !U I 

V > RAOIAN!V> 

WO = RAO|AN!WO) 

ELEMI U X A 
FLEMl 21 * E 
ELEMI 31 = I 
ELEM! 4> = OM 

“un PRECIilfING PAGE BLANK NOT FILMED 

cLFWSoJ ® WO 
X = ZERO!!! 

Y X ZER0I3I 
OX X ZER0I2) 

OY X ZER0I4) 

XO X A ♦! l.nO - E ♦ E >/! 1.00 + E * OCOSIVI ) 

LL » O 
NE » 0 
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CALL PN0L2( KUTMOL »NEO. TINCR pLL ® INTP«N6 »T I »ELEM*DELM) 
t. MRtTE(6»9» 

t 9 FORMAT* •"END OF RUN M 

99 RETURN 
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SUBROUTINE 0ERIV( T, EL, DEL ) 

IMPLICIT RFAL*8( A-H»0-2 I 
REAL*a N»NBARi,NPNfNMN»NPOf!vNSS 

REAL*R LAMl »LAM2 ,L AM3 «L AM4 « LAM5 « L AM6 , LAM7 , LAMB * LAMO vLAHlO. 

L LAM11,LAM12 

LDOICAL PRINT 

DIMENSION ELEM(200> tOELM(200) » ELI 200 ) »DEL (2 00 > vZER0(4 ) 
EQUIVALENCE lELEMI lUAI, (ELEMI 2(f E)« ( ELEM4 
F (ELEMI (ELEMt 5)»U», (ELEMI 6)»W0) 

EQUIVALENCE (OELMI 1)«DA), (DELMI 2(90E (DELM( 3 ) 901)9 
E (DELMI A)9D0M)9 (OELMI 5)«DU), (OELMI 6)90W0) 

COMMON / CONST / F J2 9 FMU 9 RSO 9 ZERO 
COMMON /FLAG/ PRINT 9 J 9 NEO 
COMMON/ ANS / X 9 V 9 XO 9 DX 9 OY 
COMMON / ALFA / ALFO 9 ALF 1 9 ALF 2 9 ALF3 

NAMELIST /TRIG/ COSV 9 S INV 9 COSU 9 S INU 9 C 0 S 1 9 SlNI 9 COS 2 U 9 SIN 2 U 9 SIN2 1 9 
U COS2I9TANU9COTI 
NAMELIST /NAM2/ V 9 U 9 I 9 P 9 X 09 N 
NAMELIST /NAM3/ DA 9 OE 9 DP 9 C 9 DV 9 WO 9 DI 9 DOM 
NAMELIST /NAM<»/ XOO 9 XOI 9 XIO 
NAMELIST /NAM5 / YOO 9 VOI 9 YIO 
NAMELIST /NAM 6 / N 9 DOCI 9 NBAR 
NAMELIST /NAMl/ 0X0090X0190X10 
NAMELIST /NAM7/ OYOO 9 OYOI 9 DY 10 

NAMELIST / lambda / L AMI , LAM2 9 LAM3 9 L AMA 9 LAM59 LAM 69 CAM 79 LAM 89 
N LAM 99 LAMIO 9 LAMII 9 LAMI 2 

NAMELIST / COEFl / A 1 9 R1 9 C 1 9 01 9 ALF 3 D 9 ALF 3N 
namelist / C0EF2 / A29B2,C2902 
namelist/ timed / TXO 9 TYO 9 TDXO 9 TDYO 

SUBROUTINE OERIV COMPUTES THE DERIVATIVES OF THE ELEMENTS ARRAY. 

INPUT 

T - TIME 

ELEM - ARRAY OF ELEMENTS 
OUTPUT 

DELM - ARRAY OF DERIVATIVES OF ELEMENTS 

DO 100 L 3 I 9 NEQ 
ELEMIL) * EL(L) 

OELMIL) » OEL(L) 

100 CONTINUE 

1 = OMODII 93.602 ) 

MO - DM0DIW093.6D2) 

IF! .NOT, PRINT ) GOTO 500 
C MRITE(692001) 

C2001 FORMAT) "0 ENTER® ) 

MRITE(69?000) ( ELEM(K)9DELMIK)9K=l9NE0) 

2000 FORMAT! 22X9 ® ELEMENT S ® 9 9X 9 ® OER ! V AT 1 VE S' // ( lOX 92020. 8 ) ) 

500 CONTINUE 

C COMPUTE FREQUENTLY USED TRIGONOMETRIC FUNCTIONS 
S!NI > OSIN! I ) 

SIN2I=nSIN(2.00* I ) 

COS I ■* OCOSII) 

COS2I=OCOS(2.00* I) 

COTIaOCOTANI 1) 

P*A*(1 ,00“E*6) 

N s OSORTI FMU / A / A / A ) 

OMO = .7500 *N «FJ2 *RSO /P /P *I 5. DO * C0SI**2 -1.00 ) 

20 V » U - WO 

XO « P / I 1.00 * E * DCOS(V) ) 

RDA « RSO /A /A 
NPC « N ♦ WO 
EPS • 1.500 ♦ FJ2 
COSV*DCOS(V) 


00004600 

00004700 


00000300 
00000400 
00000500 
00000600 
00000700 
00000800 
00000900 
00001000 
0000 non 
00001200 


00005600 


00007100 


00005600 


73 



00005900 

00006000 

00006100 

00006200 

00006300 

00006R00 


OA = -3.00 ♦ FMU * FJ2 • RSO / X0**4 /( N * OSORTt 1.00 - 6»E U 

1 ♦(( loOO -3. DO • SINI#SINI*SINU#SINU ) * E * OSIN(V) ♦ 

2 f> / XO * SINI * SINI * SIN2U ) 

DE = -EPS • FHU * RSO / X0**4 ♦ OSORTl 1.00 -E*E > / N / A *U 

1 1.00 -3.00 * SINI*SINI * SINU*SINU )♦ OSlNfVI + XO / P *( E 

2 * ncosiv) *( OCOS<V) + 2.00 > + E S1NI*SINI ♦ SIN2U ) 

01 * -EPS ♦ FHU ♦ ROA ♦ SIN2I ♦ SINU ♦ COSU / H / 

0 OSORTl 1.00 - E*E > / XO / XO / XO 
DOM s 01 ♦ TAMO / SINI 

0U= N*A*A* OSORK 1.00 - E*E ) / X0**2 - DOM * COS! 

NBAR * N ♦ DOM * COS! 
one I = DOM ♦ COS! 

MRITE(6.NAM6) 

MPN = MRAR * H 

NMN = MBAR - N 

TNPM » 2.00 » MBAR * N 

TNMN = 2.00 ♦ NBAR - N 

NSS > MBAR * NBAR * SIM! * S!NI 

OZT = ROA * NBAR *( 1,00 -1.500 ♦ SIN! ♦ S!NI ) / 8.00 
|F( PRINT ) WRITE(6,NAM3) 

A1 ® 0.500 / NBAR+ E / 2.00 ♦ NBAR / N «( N-4,00 * NBAR )/ TNPN / 

1 TNMN +3,00 / 16.00 * ROA * NSS / NBAR * EPS / 

£ (NBAR++2 - NPn**2) + EPS * nZT / 2=00 / NBAR / NBAR 

Bi = 1.00 + E ♦ NBAR / N *( A.OO «NBAR**2 ♦ N*NBAR -8. DO *N*N ) 

1 /( A.OO *NBAR*+2 - N*N ) + EPS ♦ ROA ♦ NSS / IF. DO 

2 / NPO *!-( 3.00 *NPn +5.00 *NBAR 1/«NPN+M0) 

3 /( TNPN +W0 ) + 3.00 *NP0 -5.00 *NBAR /I NPO-NBAR )/(TNPN-WOn 

Cl = 1.00 + E « NBAR *((NBAR +8.00 * N !/ TNPN +( 3.00 * NBAR 

1 -8.00 * N )/ TNMN )/ 2,00 / N + EPS ♦< OZT/ NBAR - ROA • NBAR 

2 * NSS *< 3.00 ♦ NPO + 5.00 * NBAR 1/ 16.00 / 

3 NPO /( NPN ♦ WO )**2 /( TNPN + WO ) - ROA / 8.00 ♦ NSS 

A / ( NPN + WO )**2 + ROA * NSS ♦ NBAR 

5 •( -3,00 ♦ NPO +5,00 * NBAR ) / 16,00 / NPO /( W0-NMN)**2 

6 /! TNMN - WO > - ROA * NSS / 8.00 / 

7 (WO - NMN )**2 ) 

01 = 1.00 -E ♦ NBAR++2 /N *( 9,00 * NBAR - N !/ TNMN / TNPN + EPS 

1 ♦ ROA • NSS / 2,00 / NPO / NPO ♦( NBAR * NBAR / 2.00 / ( NBAR 

2 ♦ NBAR - NPO • NPO J - 1,00 ) 

A2 = 2,00 * NBAR *( 1.00 + E ♦ NBAR / N ) - EPS * ROA * NSS / 8,00 

1 / NPO *((3.00 ♦ NPO +5.00 • NBAR ) /( TNPN + WO I +( 3.00 ♦ NPO 

2 -5,00 * NBAR )/( TNMN - WO ! ) 

B2 = E * 2.00 * NBAR ♦ NBAR * NMN / TNMN / TNPN + EPS * ROA * NSS 
1 ♦ NBAR / 2.00 / ( NBAR * NBAR - NPO * NPO ) 

C2 = -E * NBAR ♦ NMN / TNMN / TI»(PN -EPS ♦( 0.37500 « ROA * NSS / 

1 ( NBAR * NBAR - NPO * NPO ) + OZT/ NBAR ! 

02 = -2.00 * NBM *( 1.00 + E * NBAR / N > ♦ EPS •( -2.00 • OZT+ 

1 ROA * NSS *(( 3,00 * NPO +5.00 ♦ NBAR)/ ( NPN ♦ WO ) • NBAR / 

2 NPO /(TNPN + WO » - ( 3.00 ♦ NPO -5.00 ♦ NBAR )/( NMN 

3 + W0 ) ♦ NBAR / NPO / ( TNMN - WO )» / 8.00 ♦ ROA * NSS /9.00 

4 *( 1.00 /( NPN + WO )- ( WO - NMN ))) 

ALF3N « Cl ♦ ZER0(2» - A2 * ZEROIB) 

ALF30 - B2 * Cl - A2 * 01 


SINV»OSIN(V) 

CnSl)aOCOS(U) 

C0S2I)=0C0S( 2,00*0) 

SINU=OSIN(U) 

SlN2IJ=0SIN(2.00*(n 

TANtleOTANdJ) 

IF) PRINT ) WRITE(6,TRIG) 

C = enSV *) l.OO -3.00 * SIN! * SIN! * SINU * SINU )/ X0**9 
C -( 1,00 + XO / P )* SIN! * SIN! * SIN2U * SINV / X0**4 
IF( PRINT > WRITE(6,NAM2 ) 




^ ^ W IN3 »— 


c 

c 


c 

c 

c 


ALPHA TPRMS 


ALFO 

ALFl 

ALF? 

ALF3 


( 02 
( B2 
I A1 
» Cl 


ZERO! 1 » 
ZBROO » 
ZBROl^fl 
ZERO(21 


-fll 

-D1 

-C2 

-A2 


ZERO<<<tl 
ZERn<2) 
ZERfK 1) 
ZERO(3) 


>/< 

J/( 

}/l 

)/( 


Al 

fl2 

Al 

B2 


02 

Cl 

02 

r.i 


Bl 

A2 

Bl 

A2 


C2 

D1 

C2 

D1 


LAMBDA TERMS 


NBAR + N )/ 


LAMl = -< NBAR**2 +8oDO *N *NBAR )/ 2.00 *( 4,00 « 

I 7NPN * ALFl 

LAM2 = ( 3.D0 * NBAR**2 -fl.DO ♦N *NBAR )/ 2,00 ♦( 4.00 ♦ NBAR -N) 
1 / TNMN * ALFl 

LAMB = -C NBAR»*2 +8,00 •N *NBAR )/ 2,00 ♦( 4,00 * NBAR + N )/ 

1 TNPN * ALF2 


LAM4 = ( 3.00 * NBAR*'^2 -8,00 *N *NBAR 1/ 2,00 *( 4,00 * NBAR -Nl 


1 


/ TNMN ♦ ALF2 

LAMS s -NBAR / N * ( 4,00 * 

LAM6 = 2.00 * NBAR**2 / N * 

LAM7 ST ROA * NSS 
lamb = .7500 * ROA * NSS 
LAM9 = ROA * NSS 
1 / (WO + NPN 1 ♦ ALFl 

LAMIO = -ROA • NSS 
I 4,00 / ( WO - NMN ) ♦ ALFl 

LAMll = ROA ♦ NSS 

1 4,00 / ( WO + NPN ) • ALF? 


NBAR - N »/ 2,00 * ALFO 
NMN * ALF3 

* NBAR / NPO * ALF3 

/ NBAR « ALFO 

*( 3.00 * NPO +5.00 ♦ NBAR >/ 4,00 


♦( 3.00 * NPO -5.00 ♦ NBAR 1 / 
*( 3.00 ♦ NPO +5.00 * NBAR ) / 


LAM12 = ROA • NSS 
I 4.00 / ( WO - NMN ) 

IF( T.GT. 0,00 ) GOTO 300 
X = ZFROdl 
V = ZER0(3) 

OX = ZeRfl(2» 

OV = ZER0(M 
IF( .NOT. PRINT ) 

TXO = ALFO * Al + 

TYO = ALFl * Cl + 

TOXO = AlFl # A2 + 

TDVO = ALFO * C2 + 
WRITE(f.,HMEO) 

GOTO 400 
300 CONTINUE 

XOO ® ALFl * DSIN( 

I + ALFO / 2.00 


♦( 3.00 ♦ NPO -5.00 * NBAR ) / 


* ALF2 


GOTO 400 
ALF2 • Bl 
ALF3 • 01 
ALF3 * B2 
ALF2 * 02 


(X:oS( 2.00 ♦ NBAR * T > 


XOl = LAM6 / TNMN 


1 


XIO 


2,D0*NflAR#T » +ALF2 * 

/ NBAR 

/ TNPN ♦ OS INI N*T ) + LAMS / TNMN / TNPN • 

DCnS( N*T I - LAMl / N /( 4,00 4NBAR +N 1 ♦ DSINI TNPN*T I + 
LAM2 / N /( 4.00 *NBAR -N 1 * 0SIN( TNMN*T 1 -LAM3 / N / 

I 4,00 *N6AR + N 1 4 OCnSI TNPN*T I + LAM4 / N / ( 4,00 *NBAR 
-Nl+OCOSI TNMN *T ) 

= ( LAM7 * DSINI 2.00 ♦ NPO • T ) ♦ LAMB * OCOSI 2,00 * NPO * 

T n / 4,00 / I NBAR * NBAR - NPO ♦ NPO ) - I LAM9 ♦ DSINI 


2,00 *( NPN + WO »* 
n / 4,00 / NPO / I 


T ) + LAMll • 
TNPN + WO } + 


OCOSI 2,00 *( NPN ♦ WO 1* T 
I LAMIO * OSINI 2,00 *I WO 


YOO 

VIO 


- NMN >♦71+ LAM12 • OCOSI 2.00 ♦( WO - NMN 1* T >>/ 4,00 
/ NPO /I TNMN - WO 1 + ,500 « OZT * ALFO / NBAR / NBAR 

* ALFl ♦ OCnSt 2,00 * NBAR ♦ T ) - ALF2 ♦ DSINI 2,00 ♦ NBAR • 

T > ♦ ALF3 

t OZT / NBAR *( ALFl * OCOSI 2,00 * NBAR * T 1 - ALF2 ♦ OSIN 

I 2,00 ♦ NBAR « T n + NBAR / 4,00 / NPO /( NBAR ♦ NBAR - 

NPO * NPO i *( LAM7 ♦ OCOSI 2,00 * NPO ♦ T > - LAMB * DSINI 
2,00 * NPO * T >1 - RDA * NSS ♦ ALF3 / 2,00 / NPO / NPO • 

OCOSI 2,00 * NPO * T > + NBAR / NPO / I NPN + W0 ) / 4,00 / 

I TNPN + WO t»I LAMll * DSINI 2.00 *1 NPN + W01«T)- LAM9 « 
OCOSI 2,00 *1 NPN + WO m 4 RDA * NSS / 8.00 ♦IIALF2 * OSIN 
I 2,00 *( NPN + WO )*T> - ALFI * OCOSI 2,00 ♦( NPN + WO )* T 

n/( NPN ♦ HO >**2 - I ALFI * OCOSI 2.D0 • 
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9 C WO - NMN )♦ T J + ALF2 * DSIN( 2„0O ♦( WO - NMN >♦ T >) / 

A f WO - NMW )**2 ) + NBAR t NPO / <».DO /( WO - NHN ) /( TNMN 

B - WO lAMlO * OCOS( 2*00 *( WO - NMN )* T ) - LAM12 • OSIN 

C I 2. DO *( WO - NMN I* T )) - 02 T * ALFO * T / NBAR 

YOl = - NBAR / N •( ALFO / N * OSINI N*T I + NBAR / N • ALF3 • 

1 DCnSI N*T I + 2*00 / TNPN / TNMN *( LAM5 * OSINI N*T ) - LAM6 

2 * OCOSI N*T >))+ I *>*00 ♦ NBAR / N /( TNPN *2*00 • NBAR I ♦ 

3 ( LAM3 * DS!N( TNPN*T ) - LAMl * OCOSI TNPN*T M - NBAR * 

A ( NBAR +8*D0 * N )/ TNPN •( ALF2 • DSINf TNPN*T I - ALFl * 

5 DCnSI TNPN*T tl) / 2*00 / TNPN +( NBAR ♦( 3*00 ♦ NBAR -B*DO 

6 * N )/ TNMN ♦( ALF2 * OSINI TNMN*T I - ALFl ♦ DCOSI TNMN*T H 

7 ~ A. DO NBAR / N / I TNMN +2,00 * NBAR)* ( LAMA * OSINI TNMN 

8 *T ) - LAM2 • DCOSI TNMN*T >11/ TNMN **500 

DXOO = 2*00 * NBAR *IALFI ♦ DCOSI 2.00 * NBAR ♦ T I - ALF2 * 

1 OSINI 2.00 « NBAR * T I) 

0X01 a N / TNMN / TNPN *I LAM6 * DCOSI N»T I ~ LAMS ♦ OSINI N*T I) 

1 - TNPN / N /I A, DO * NBAR + N 1*1 DCOSI TNPN ♦ T I ♦ LAMl - 

2 OSINI TNPN • T )) ♦ LAM3 + TNMN / N /I A*00 • NBAR - N I* 

3 I LAM2 * DCOSI TNMN • T ) - LAMA ♦ OSINI TNMN * T I) 

0X10 = 0,500 • NPO /I NBAR*NBAR - NPO*NPO 1*1 LAM7 ♦ DCOSI 2. DO * 

1 NPO *J I -LAM8 • OSINI 2*00 ♦ NPO * T M+ 0.500 *1 WO + NPNI 

2 / NPO /I WO + TNPN I*I LAMll * 0SINI2.D0*! W0+ NPN 1*T I 

3 - LAM9 * DCOSI 2,00 *1 WO+NPN >*T I) +0*500 *( WO - NMN I / 

A I TNMN - WO I^l LAMlO ♦ DCOSI 2*00 ♦! WO«NMN l*T I - LAM12 

5 • DSINI 2,00 WO-NMN >*T 11/ NPO 

OVOO = -2*00 * NBAR *1 ALFl ♦ OSINI 2.00 ♦ NBAR * T I ♦ ALF2 * 

I DCOSI 2*00 * NBAR ♦ T H 

OYOl a -NBAR *I ALFO / N +2.00 * LAM5 / TNMN / TNPN I* OCOSI N*T I 

1 +N8AR A I NBAR / N * ALF3 -2.00 ♦ UAM6 / TNPN / TNMN ) * 

2 OSINI NAT I + NBAR « 0.500 AI A.OO * LAM3 / N /I A. 00 * NBAR 

3 + N ) -I NBAR +8,00 + N 1/ TNPN * ALF2 I* DCOSI TNPN ♦ T I 

A + NBAR * ,5D0 ♦!! 3,00 * NBAR -B.DO ♦ N )/ TNMN * ALF2 -A.OO 

5 ♦ LAMA / N /I A.OO * NBAR - N ))♦ OCOSI TNMN*T » + NBAR ♦ 

6 0*500 *I A.OO * LAMl / N /I A.OO » NBAR + N ) -I NBAR +8.00 

7 AN)/ TNPN A ALFl JA DSINI TNPN*T I + NBAR A 0.500 *11 3.00 

8 A NBAR -RoDO * N )* ALFl / TNMN -A.OO * LAM2 / N /I A.OO * 

9 NBAR - N M* DSINI TNMNAT ) 

OYIO = -2,00 A 02T A| ALF2 * DCnSI 2*00 * NBAR * T ) + ALFl * 

1 DSINI 2. DO A NBAR * T II - NBAR * LAMB A 0,500 /( NBARANBAR 

2 - NPOANPO )A DCOSI 2,00 • NPO * T ) -( NBAR * LAM7 * 0,500 

3 /I NBARANBAR - NPO*NPn ) - ROA * NSS * ALF3 / NPO ) ♦ OSINI 

A 2, DO A NPO A T )+I LAMll * NBAR * 0,500 / NPO / I TNPN + WO) 

5 +0,2500 A ROA * NSS * aLF2 /( NPN + WOII* OCOSI 2.00 *1 NPN 

6 + WOI* T ) -I LAM12 A MBAR * 0,500 / NPO /( TNMN - WO I ♦ 

7 ,2500 A rdA a NSS • ALF2 /!• WO - NMN ))* OC0SI2,00 ♦( WO-NMN 

8 )AT I +1 LAM9 A NBAR AO, 500 / NPO /I TNPN + WO >+ 0.2500 * 

9 RDA A NSS A ALFl /I NPN + WOIIA DSINI 2,00 * I NPN+WO) *T )-l 

A LAMlO A NBAR aq.SDO / NPO /I TNMN-WO I -0.2500 * RDA A NSS 

B A ALFl /IWO - NMN )IA 0SIN(2,00 *I WO-NMNIATI-OZTAALFO / NBAR 

X * XOO + € A XOI + EPS A xlO 
Y = YOO + E A YOl + EPS * VIO 
OX » 0X00 + E A 0X01 + EPS * 0X10 
OV « OYOO + E A DYOl + EPS * OYIO 
400 IFI .NOT .PRINT I GOTO 600 
WRITEI6,C0EFU 
WRITEI6,C0EP2t 
WRITE (6,NAM4I 
WRtTEI6fNAM5I 
WRITEI6.NAM1 1 
WRITEI6,NAM7I 
WRITEI6, LAMBDA! 

600 CONTINUE 

00 200 L>lcNEO 
DELU) a OELMIL) 

200 CONTINUE 

99 RETURN 00012500 

END 
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SUBRnuTINE OUT) TvELEMtOjERROR »N»LK«M ) 

implicit real*bu-h«o-z > 

REAL*® I 
LOGICAL PRINT 

DIMENSION T(1 ) »ELEM< n ,D( 1 ) «ERRDRI 1) 

DIMENSION DEG(5)«ELE(S00I 
COMMON /FLAG/ PR!NT«J,NEO 
common/ ANS / X*V»XO»DX(DY 
COMMON / ALFA / ALFO, ALF 1 , ALF 2 » ALF 3 

EOIMVALENCF (ELE I U«AU (ELE ( 2I«EK (ELE ( 3)«n. 

E (ELE ( 4I»OMl, (ELE ( SKUl, ( El.E ( 6KW0) 

DATA TPI/ 6.2031853072 / 

OEGREE(X) = X • 57.2957795130823 
VMAG(X,y» « DSORT( X*X + Y*Y ) 

J * J ♦ 1 
URITE(6«3> 

3 FORMAT(130(lH-» ) 

POS => VMAG(X,Y f 
DO 10 Kal»NEO 
ELE(K) = ELEM(KI 

10 CONTINUE 

VEL » VMAG(DX»OV) 

WRITF(6»4I POSoVEL 

4 FORMAT! »0S25X* 'POSITION = S 020. 10, 15 X, • VE LOC I TV = ',020.10 ) 

V = U - WO 
DO 20 K=3,6 

20 ELE (K) = OMOD( FLE !K»,TPI ) 

V = DMOD( V,TPI I 

0EG(1» = DEGREE! ELE (3> ) 

DEG(2I = DEGREE! ELE (A)> 

0EG(3) = DEGREE! ELE !6I i 

0EG(«) = DEGREE! ELE (5M 

DEG(5) - DEGREE! V ) 

WRITE(6,1001) T,X,Y,XO 

1001 FORMAT! ,4X, 'TIME (SECI = • ,015.8, ®X, «X !KM) = ',O15.0,8X, 

1 *Y !KM) = 9 ,D15.8,8X, 'XO (KM) = ',015.8 ) 

WRITE!®, 1003) A,E,1,DEG!1) 

1003 FORMAT! 'O', ex, 'A (KM) = ' , 015 .8 ,9X , ' ECCEN = ' ,015.8, lOX, ' INCL « * 

3 ,015.8,9 RAO - 9,015.8,9 OEG') 

WRITE (6, 1004) ELEM!4) ,DEG!2) ,ELEM!6) ,DEG!3) 

1004 FORMAT! «0» ,8X , «L .A. N. = ',015.8,9 raq n 9,015.8,' DEG9,10X, 

4 'ARG OF PER * ',015.8,9 RAO » 9,015.8,9 OEG') 

WRITE!®, 10051 V,DEG(5 I ,U,0EG!4) 

1005 FORMAT! 'O', 13X,«V = ',015. 8,' RAD = ',015.0,' 0EG',19X,'U = ', 

5 015.8,9 RAO = 9,015.8,' OEG') 

WRITE !6,1006) ALFO, ALF 1 , ALF2 , ALF3 

1006 FORMAT! 'O' ,10X, 9 ALFO = « 015 .8, I IX , *ALFl ='015.8, IIX, «ALF2 >'015.8, 
i 11X,'ALF3 =*015,8) 

50 RETURN 
END 
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SUBRnilTINE TERM( T,Et.EM,DFL 

IMPLICIT RFAL*fl (A-H,n -2 ♦») 

dimension fLEMUI* DELMdl 

COMMON /INTPG/ TFoTINCRtTd 

COND e TE - T 

RETURN 

END 


MfCOND » 
INTP»KUTM0L 





BLOCK OATA 

IMPLICFT REAL*8( A-H,0-2*S F 
RFAL*8 f 
LOGICAL PRINT 
OIHFNSION 2FR0(4> 

COMMON / CONST / FJ2,FMU.RSO*ZERO 
COMMON /INTFG/ TF » STFP , T I , IN TP ,KUTMOL 
COMMON / ELMFNT / A ^E I vOMtV 9U«W0 
COMMON /flag/ PRINT, JvNEO 
OATA FJ2 / 1,08230-3 / 

OATA FMU,RSQ / 3.9S6032D5, 4o06fl09fl07/ 

OATA ZFRO / 4*1,00 / 

DATA STEP,TI ,TF,rNTP,KUTMOL/ 6.02,0,00,8,6404,1,2 / 
DATA A,E,I,0,V / 1.500, 1,0-2, 45,00, 2*0,000 / 

DATA OM,WO / 2*45.00 / 

DATA PRINT /.FALSE./ 

END 


li 

Vi 




//gObOata's nn ♦ 
filNIT 
Fj?*O.DO, 

Aa?cODA» 

F = 1 .D- 4 » 

t “10 oDO(t 

nM“ 20 . 00 <i 

V ®45 oOOe 
WC“0.00» 

2F^n«l ,D0,1.,D-3,1 o'OOtloD-3 

CFwn 

filNTEG 

TF*fe,04Bn5» 

NE 0“69 

INTP=6« 

CFNn 


